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PREFACE 


Our interest in scattering by nonspherical particles was originally motivated by a desire to understand 
the impact of atmospheric aerosols on the Earth’s climate. Since such aerosols are routinely assumed 
to be spheres, yet few except haze drops are spheres, we were naturally interested in what the dif- 
ferences in scattering and absorption might be. 

After a time, however, we realized that our study had a considerably wider applicability. Scattering by 
nonspherical particles is important in practically every field of science and engineering. Biologists 
study nonspherical cellular structures by light scattering. Medical laboratories analyze blood cor- 
puscles. Astronomers study interstellar grains. Engineers analyze coal slurries. And on and on. What 
all these enterprises share, is a deep lack of knowledge of nonspherical scattering. Hence we began to 
see our calculations as of much wider significance. Specialists in all these fields may find them of 
value. 
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I. INTRODUCTION 


When we first embarked on a course of research into the scattering properties of nonspherical par- 
ticles, our goal was simple: to search for the main differences between spherical and nonsperical parti- 
cle scattering by examining a number of special cases. By taking such an inductive approach, we essen- 
tially cast ourselves in the role of experimenters, with the important distinction that we solved Max- 
well’s Equations on the computer rather than letting Nature do so in the laboratory. This allowed us 
to obtain data which is very difficult to measure, and hence rarely measured — for example, the ab- 
sorption cross-section and the backscattering near 180 degrees. 

Our initial motivation was not to examine the scattering features of individual particles, but rather to 
unearth scattering features common to all nonspherical particles. In order to do so, we selected what 
we considered to be a rather general class of particles which could be continuously deformed from a 
■sphere. These “Chebyshev particles” are obtained by rotating the curve 

r s - r 0 [1 + e T n (cos 0 S )] 

= r 0 [l + e cos(n© s )] (I) 

about the vertical axis 0 S = 0 (T n is then the n-th Chebyshev polynomial and e is the deformation 
parameter). Various examples of these particles are shown in Figure 1. 

Following an extensive literature survey of scattering methods, summarized in Appendix A, we 
selected the “Extended Boundary Condition Method” (EBCM) invented by Waterman (1965, 1971) 
and further developed by Barber and Yeh (1975). Prof. Barber of the University of Utah was kind 
enough to supply us with a copy of his EBCM computer code. After considerable modification to im- 
prove its efficiency and speed, and to automate the convergence testing, we ran a large number of 
cases: 23 different particles in all, with deformation parameters e = -0.20 to 0.20 (in steps of 0 . 05 ) 
and ‘equal-volume-sphere’ size parameters x = 1 to 25 (in steps of 1). 

Our initial paper (Mugnai and Wiscombe, 1980) examined only a tiny subset of all this data — just the 
extinction and absorption cross-sections and the phase function at 180 degrees for size parameters x 
< 10. We reserved the study of phase function and degree of polarization to a future paper. As we 
planned this second paper, however, it became clear that we would still be able to show only the 
smallest fraction of our data, namely that part which was most indicative of general trends. This, of 
course, was in line with our original goal. 

But the question naturally arose, should we throw out the remainder of our data? Here the analogy 
between our study and an experimental investigation pressed itself upon us. A good experimenter will 
publish some of his basic data as well as his interpretation of it. The reason is, that others may be able 
to use the data, or draw conclusions from it, in ways not envisioned by the experimenter. 
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ly, these selection criteria embodied some preliminary data interpretation, but no further analysis of 
the results is undertaken here; that is reserved for a forthcoming journal article. 


Appendix A, which is condensed from an intended review article, summarizes the entire field of 
nocspherical scattering. The book edited by Schuerman (1980) and our first paper (Mugnai and 
Wiscombe, 1980) can also be recommended as surveys of the state-of-the-art. 



II. PARTICLE SHAPE 


It seems to us that, along with ellipsoids and cylinders, Chebyshev particles constitute the most in- 
teresting class of fundamental nonspherical shapes whose scattering properties deserve attention. 
Among all simple analytical forms which we considered, they seemed the most attractive for several 
reasons: 


• they differ distinctly from other shapes for which extensive calculations have been made 

(spheroid and cylinder) 

® by increasing the magnitude of e, a continuous deformation away from a sphere can be 
achieved — a property Chebyshev particles shared with spheroids, but not, for example, with 
cylinders and cubes 

• while starting as convex, they quickly develop concavity as the magnitude of e increases 

• they exhibit surface roughness in its simplest possible manifestation, namely “waves” of 
uniform amplitude running completely around the particle; the effect of surface roughness 
can be studied in a systematic way, by varying e and the waviness parameter V 

8 the surface of the particle is smooth in the mathematical sense — all its derivatives exist and 
are bounded; hence highly-efficient Gauss quadrature rules can be employed at key points in 
the EBCM calculation 

• the particles are rotationally symmetric, the only practical case for the EBCM (see Barber, 
1973, Dp. 43-50; also Mugnai and Wiscombe, 1980) 

• The Chebyshev polynomials form a complete set, so that any rotationally-symmetric shape 
can be represented as a sum of elementary shapes of the form of Eq. (1). 

Although we were unaware of their work when we selected Chebyshev particles, Pruppacher and Fit- 
ter (1971) assumed the shape of deformed raindrops to be a sum of elementary Chebyshev particles. 

This is a natural choice, being none other than a Fourier cosine series. 

The 23 Chebyshev particles which were used in this investigation are: 
for e = ±0.05: T 3 , T 4 , Tg, Tg, 120 
for e = ±0.10: T 2 , T 3 , T 4 , T 6 , T g 
for e = ±0.15: T 3 , T 4 
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for e = ±0.20: T 2 

Tj particles are the same whatever the sign of e, hence there are only 23 different particles in this seem- 
ing list of 26. A 24th ‘particle’, which we call MIXTURE, is a blend of these 23 in equal proportions. 
Figure i presents 3-D pictures of all these particles except the ones with e = ± 0.05 (which do not look 
much different from a sphere). 

We shall refer to these particles with the following notations, the first of which was used in Mugnai 

and Wiscombe (1980): 

T n + , T n ~: Chebyshev particle in fixed orientation with waviness parameter V and 

positive ( + ) or negative (-) value of deformation parameter e 

T n ( + e), T n (-e): Chebyshev particle in fixed orientation with waviness parameter ‘n’ and 
deformation parameter e = ± e or - e 

For Chebyshev particles in random orientation, we use a similar notation, but now between angle 

brackets, e.g. <T n ( + e)>. 

T 2 (±0.10) particles are convex and quasi-spheroidal (prolate for +, oblate for -). T 2 (±0.20) par- 
ticles, however, have almost perfectly flat sides ( + ) or perfectly flat tops and bases (-), more 
resembling cylinders or disks (with rounded edges) than spheroids. For even larger deformations, the 
T 2 particles become concave, with T 2 + tending toward dumbbells and T 2 ~ toward toroids. 

T 3 particles are pear shaped, T 3 + and T 3 ~ are actually the same particles; this is true for all T n par- 
ticles with odd ‘n’ (but not with even ‘n’). They are the only case we study that is not mirror-symmetric 
in the plane perpendicular to the rotation axis (there are substantial computational savings in the 
EBCM, and different branches in the code, for particles with such a plane of symmetry). To our 
knowledge, ours are the first published EBCM results for particles lacking a plane of symmetry since 
Waterman (1971) studied a cone-sphere combination. 

T 4 + particles have two shallow valleys girdling them at mid-latitudes; T 4 ~ particles, one girdling valley 
at the equator, plus two dimples at the top and base. T 6 and T g particles are similar, but with more 
valleys; the only new features are pimples (instead of dimples) at the top and base for the ‘ + ’ cases. 
All these particles have either bulges at the waist, or pinched waists. 

In spite of their notational similarity, T n ( + e) and T n (-e) particles look quite different to the eye. 
They bear a subtle ‘conjugacy’ to each other, however; for if we imagine the sphere r s = r 0 as a mir- 
ror, the two particles are mirror reflections of each other. It is therefore of some interest for the reader 
to compare scattering quantities for these particle pairs. For that reason, their plots are made nearest 
neighbors in the various plot groupings to follow. 



PARTICLE SHAPE 


The particle refractive index was fixed at m = 1 .5 - 0.02 i, a value typical of some maritime aerosols 
in the visible region (Gerber, 1979). In Mugnai and Wiscombe (1980), we explored imaginary indices 
of 0, 0.05 and 1.0, but for the present study, we felt it was more important to explore as large a range 
of particle shapes and sizes as possible. We avoided a zero imaginary index because, in that case, the 
Mie results are full of the worst spikes and “ripple” (see Sec. VIII). This jagged behavior of the Mie 
quantities is reflected, in somewhat muted form, in the nonspherical quantities, making it more dif- 
ficult to obtain convergence in the EBCM, with less accurate results. On the other hand, for very large 
imaginary indices, the spherical-nonspherical differences become very small and therefore less in- 
teresting, even though EBCM convergence is excellent. 

All our EBCM results for Chebyshev particles are compared to results for an equal-volume sphere. 
Let V n (e) be the volume of a Chebyshev particle as given by Mugnai and Wiscombe (1980) (note that 
their expression for particle surface area is incorrect; the correct formula is given by Chylek, et. al., 
1982). Then the radius of the equal- volume sphere, r ev , is defined by 

V n (e) = 4 7r r ev 3 /3 ( 2 ) 


and the equal-volume-sphere size parameter, x, is given by 

x = 2 7 r r ev /A ( 3 ) 

where X = wavelength. The units of r ev and X are neither needed nor specified in this study; only the 
non-dimensional parameter x is used to characterize the size of the scatterer. 

In our calculations, we always chose a value of x a priori, then adjusted the r 0 of each Chebyshev par- 
ticle in order to yield that value of x (actually, it is r 0 /A that is adjusted). Thus, for each fixed x, the 
‘basal radii’ r 0 of the various Chebyshev particles differ slightly. In particular, r 0 will be different for 
T n + and T n ”" particles having the same size parameter x. Thus, this pair of particles will not be exactly 
‘conjugate’ in the sense discussed above. 




T 2 (-0.10) 


T 2 {-0.20) 


Figure 1. Three-dimensional drawings of 18 of the Chebyshev particies used in this 

scattering study. 














Let us further define an (arbitrary) reference axis z’ through the scatterer. (For axisymmetric particles 
such as ours, the problem simplifies immensely if the reference axis and the rotation sods are coinci- 
dent.) Then the orientation of the particle in the laboratory frame is given by the zenith angle 
0 p and the azimuth angle <f> p of £ (see Fig. 2a again). 






EBCM — THEORETICAL CONSIDERATIONS 


E s (r) 


E D„ fe M 3 „ (kr) + b* N 3 „ (kr)] (4b) 

v = i L J 


where k = 2 tt/X; r is the radius vector; v is a combined index incorporating the usual spherical har- 
monics indices a, m,n (Morse and Feshbach, 1953, p. 1865), where a is either ‘e’ (even) or V (odd); D, 
is a normalization constant 


(2n + 1) (n-m)! 

” m n(n+l)(n + m)! 


( 5 ) 


(e 0 = 1, e m = 2 otherwise); a„' and b^ are the known expansion coefficients of the incident field; a„ s 
and b„ s are the unknown coefficients of the scattered field (to be determined by the EBCM); and M s N 
are the vector spherical harmonics 


M„ (r) = Vx 


I cos(m <j>)) 
\sin(m 4>)f 


P n m (cos 0) z n (kr) 


N„ (r)=ivxM„ 


(6a) 


(Morse and Feshbach, 1953, p. 1865). Here r, 0, and <t> are the usual spherical coordinates of the field 
point r; cos(m <j>) is used when the index a is ‘e’ (even), sin(m <f>) when it is ‘o’ (odd); P n m (cos 0) is the 
associated Legendre function; and z n (kr) is a spherical Bessel function. 

M 1 and M 1 must be finite at the origin r = 0, and therefore 


z,, (kr) = j n (kr) (7a) 

where j n is the spherical Bessel function of the first kind. M 3 and N 3 must represent outgoing waves at 
infinity, and therefore 


z n (kr) = h n W (kr) = j n (kr) + i y n (kr) (7b) 

where h n W is the spherical Hankel function of the first kind, and y n is the spherical Neumann 
function. 

We assume that M 3 and N 3 may be used to represent the scattered field even for concave p articles, for 
which the scattered wave is not entirely an outgoing wave near the particle; by so doing, we apply the 
“Rayleigh hypothesis” (Millar, 1969, 1973), which had its origin in the theory of rough reflecting sur- 
faces. 
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Since the incident field is assumed to be a plane wave, we have 

E‘ = E° e ikz (8) 

E° is a unit vector along the direction of polarization (i.e., we normalize the incident field to unity). 

Theoretical solution 


A thorough formulation of the EBCM for dielectric scatterers can be found in Barber and Yeh (1975). 
Basically, their method of solution makes use of Schelkunoff’s (1943) equivalence theorem, which 
says that the field scattered by a particle can be exactly reproduced by ‘equivalent’ (albeit fictitious) 
electric and magnetic currents on its surface. Then, applying the usual electromagnetic boundary con- 
ditions (continuity of the tangential components of the fields), and truncating the expansions (4) after 
N terms, one obtains a system of linear equations relating the scattered and incident field coefficients: 



(9) 


This is a 2N x 2N system. The expansion coefficients of the scattered field are obtained by solving it. 
In Eq. (9), the so-called “transition matrix” is given by 


K’ + rhJ’ 

L’ + ml’ 


K + fhJ 

L + ml 

I’ + mL’ 

J’ + ifiK’ 


I + mL 

I + mK 


( 10 ) 


where m is the refractive index of the particle. I, J, K, L and I’, J’, K’, L’ are two-dimensional in- 
tegrals over the particle’s surface, which define the particle from an electromagnetic point of view. For 
example, the surface integrals I are given by 

I =—■ J s 1VP (kr s ) x M* (mkr s ) • dS (11) 

Here, r s is the radius vector from the origin to the surface element dS; and the combined index v’ in- 
corporates cr’,m’,n’. 

The surface integrals J, K, and L have a form similar to (1 1), but with different vector cross-products: 
M 3 x N 1 , N 3 x M 1 , and N 3 x N 1 , respectively. The surface integrals I’, J’, K’, and L’ have the same 
form as I, J, K, and L, respectively, except that all the vector spherical harmonics are of the first kind. 

It is advantageous to calculate the surface integrals in the body frame because, in that frame, they are 
independent of the direction of the incident field. Thus, when we solve the system (9) in the body 
frame, we avoid recalculation of the transition matrix T when the orientation of the particle is chang- 
ed. Moreover, in the next sub-section, we will see that for axisymmetric particles (like ours) the com- 
putation of surface integrals in the body frame can be dramatically simplified. 


EBCM — THEORETICAL CONSIDERATIONS 
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Surface integrals in the body frame 

As an example of the formulas for the surface integrals in the body frame, let us consider: 


where 


JemnomV d ^ S ] K W ( kr S ) j ( ffikr s> ( kr s) 2 

n O _j 

x [ m Amnm'n’ ( cos © s ) sin(m<£ s ) sin 
+ (cos 0 S ) cos(m <£ s ) cos (m> s ) 

d(cos 0 S ) 


An’n’mn (COS 0) = ^ ^ Pg (COS 0) 


n -m + 1 * m , rw 

— P . (COS0) 

n +1 n+i v ’ 


n +m 


P m ,(cos 0) 
n n-i 


( 12 ) 


(13) 


and 


P n m (cos 0) 


P n m (cos 0) 
sin 0 


(14) 


Here, P n m is the associated Legendre polynomial. The other surface integrals have similar expressions, 
which can be easily obtained from Eq. (11) and the comments following it. 

For axisymmetric particles, it is possible to perform the integral over azimuth <j> s analytically, since r s , 
the radius of the particle surface, is a function of cos 0 S only (e.q., Eq. 1 for Chebyshev particles). As a 
consequence, many of the surface integrals vanish, depending on the relative parity of the indices a 
and a’, and on the values of the azimuthal indices m and m\ For example, for the surface integrals I, 
we have (Barber, 1973, pp. 43-45): 


I, 


omnom’n 


, = 0 


for all m and m 5 ; and 


*emnom n 


*omnem n 


for m’ t m. 
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Therefore, only the surface integrals I emnomn . and I omnemn , are non- vanishing. And because of the fur- 
ther simplification that 


I 


omnemn 5 


*emnomn 


(16a) 


we need only calculate 

Iemnomn’ “ m J h n W(kr s ) j n , (mkr s ) (kr s ) 2 (16b) 

-i 

Finally, for axisymmetric particles with a plane of symmetry perpendicular to the axis of rotation (like 
Eq. 1 for even ‘n’), by virtue of the following property of associated Legendre functions: 

P n m (-cos 0) = (— i) n " m P n in (cos 0) (17) 

we have 


s> + Amn’mn ( COS ©s>] d (c° s ©s) 


if n and n* are both even, or both odd; otherwise, I emnomn . is still given by (16b), but the integral can be 
done just over the half-range (0,1) and doubled, if desired. 

For axisymmetric particles, analogous simplifications exist for the other surface integrals as well 
(Barber, 1973, pp. 45-46). 

Evaluation of the scattered field in the body frame 

In order to evaluate the scattered field in the body frame, one must first express the expansion coeffi- 
cients of the incident field in that frame. 

The incident field coefficients are given by (Barber, 1973, pp. 137-139): 

a Lm = in Vn(in +T) E° • C 0 . mn (0,. 0,) (19a) 

b Lnn = -in+ 1 v»(n+i) E° • (©J, 0.) (19b) 

Here, E° is the unit vector in the direction of polarization of the incident electric field (8); 0; and 0, are 
the zenith and azimuth angles of the direction of propagation of the incident field (in an arbitrary 
coordinate system); and the formulas for the vector functions C and B are given by Morse and 
Feshbach (1953, p. 1898). 
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Let us work in the body frame (0; = 0 p , = 0), and define the auxiliary function 

f mn (x) = n x P n m (x) - (n + m) P m n _j (x) (20a) 

Then the incident field coefficients become 

a* = -i n f (cos 0 ) E° , 

emn mn v p' y’ 

bi emn = + ' f mn C COS 0 P> E V > 

ai omn= i n m P n m (cos 0 p ) E° x , , 
b* n = -i n + 1 m P™ (cos 0 J E° , 

omn n v P 7 y 

where the components of E° along the x” and y’ axes are: 

E° x> , = - E° x cos <t> p - E° y sin <f> p 

E° y , = E° x sin <j > p - E° y cos <j> p 

and E° x and E° y are the components of E° in the laboratory frame. 

Since we know the surface integrals I, J, etc., and therefore the transition matrix, we can now merely 
apply the linear transformation (9) to obtain the expansion coefficients a s and b s of the scattered field 
in the body frame. 

The vector far-field (kr —■<») amplitude of the scattered field F is defined by 

E s (kr) - F (e s ,e z ) e ikr /kr (22) 

where e z , e s are unit vectors in the directions of the incident and scattered waves, respectively. The vec- 
tor spherical harmonic expansion of F in the body frame is given by Barber (1973, pp. 148-150). Ab- 
sorbing a factor [n(n + )}' A into the spherical harmonics, the components of F parallel (q = 1) and 
perpendicular (q = r) to the plane defined by z’ and s are 

N 

F b ^ = E r (n + *> D, [a* q (0 b , <A b ) + i b) (0 b , 0 b )] (23) 

>/ = i 


(20b) 

(20c) 

(20d) 

(20e) 

(21a) 

(21b) 


© b and <j> b , the zenith and azimuth angles of the direction of scattering e s in the body frame, are shown 
in Fig. 2d and are given by: 
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(9) = tarr 1 


([(sin 0 cos 0 p cos <£ p -cos 0 sin 9 p ) 2 + (sin 9 sin 0 p ) 2 ] 1/2 1 
sin 0 sin 0 p cos <j> p + cos 0 cos 0 p 

{ sin 0 sin <£_ ) 


(0) = tan -1 


I cos 0 sin 0 p - sin 0 cos 0 p cos <j> p 


(24a) 


(24b) 


is the scattering angle in the laboratory frame). The modified vector spherical harmonics are 


C inn = B'emn = </> b ) P n m (COS© b ) 

(25a) 

C in = B r omn = m COS(m<(> b ) P n m (C0S© b ) 

(25b) 

C r emn = -Binn = -cos(m<p b ) f mn (cos0 b ) 

(25c) 

Cr omn = -B'tmn = f mn (cOS© b ) 

(25d) 


The scattered field in the laboratory frame 

The final step is to calculate the two components of F in the laboratory frame — one (F^) parallel to 
the scattering plane, the other (F r ) perpendicular .to this plane: 


F^(0) = F b X (0 b , 0 b ) cos 0 + F b r (0 b > 0 b ) sin 0 
F r (0) = -F b ^(0 b . <£ b ) sin (3 + F b r (0 b , 4> b ) cos 0 


(26a) 

(26b) 


where 

0 = tarn 1 | (cos© cos9 p cos <j> p sin 4> b + sin0 sin0 p sin<£ b (27) 

+ cos© sin <f> p cos</> b )/(cos0 p sin<£ p sin<£ b - cos </> p cos<t> b ) | 

The two components of the intensity of the scattered field (using the same notation as in Mie theory) 

are then simply 

i q (9)= |F q (0)| 2 (28) 

where q = i or r. The scattering and extinction cross-sections are given by (Waterman, 1965): 
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*ext =T?Im [E°* F(0)] 


where F(0) is the far-field amplitude of the scattered electric field in the forward direction (6 = 




For axisymmetric scatterers, all the surface integrals vanish for m = mfr (recall that m and m ? are com- 
ponents of the tripartite indices v and v’ of the scattered and incident field expansions). As a conse- 
quence, by re-ordering v = (<r,m,n) to v = (a, n,m) (and similarly for */’), and allowing a and n to run. 
through their entire sequence before m is incremented, the transition matrix T in (9) becomes block 
diagonal. Therefore, due to the linearity of the equations, (9) splits into (n max + 1) independent sub- 
systems, each corresponding to a different azimuthal mode m (m = 0. . .,ii max , where n max is the 
highest n- value used); 



Here, T m is the appropriate block of T for the given m-value, and the expansion coefficients a and b 
are for all possible combinations of a and n (or o' and if), holding m (or nf = m) constant. In what 

two indices a and a (or a ’ and n ! ). 

Because of the linearity of the vector spherical harmonic expansions, the splitting over separate 
m-values carries through the entire calculation; thus, the entire scattering problem can he replaced by 
(n ma x + 1) independent sub-problems, each corresponding to a different azimuthal mode m. This is 
accompanied by a re-ordering of the sums over m and n in all spherical harmonic expansions, as 
follows: 



n as o m — 0 


n max 




(31) 


This means that the total scattered 


can 


calcu! 


re sum 


partial fields 


fs = V 

wand 


m = 0 


(32) 


where each partial field is obtained from an analogue of (4b), but with coefficients obtained from (30), 
and a summation over a rather than v. 

This result brings important simplifications and speed enhancements in the numerical solution. For, in 
the general case, the linear equations to be solved form a potentially very large 2Nx2M system (Eq. 9), 
where 






; A: of 
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N = n max (n max + 3) (33) 

is the total number of terms to be summed in the spherical harmonic expansions (4) (note that the n = 0 
terms vanish and therefore are not taken into account when computing N). However, for axisyrn- 
metric scatferers, this system is replaced by the (n max + 1) sub-systems (30), each of dimension 2(n max 
- m + 1). Note that the sub-system dimensions decrease as m increases, because, in the matrices T m , n 
ranges from m to n max . Also, the sub-system dimensions already account for parity relations among 
the surface integrals (e.g., 16a), which reduce the size of the matrices T m by a factor of 4. 



V. EBCM — NUMERICAL CONSIDERATIONS 


There are four major numerical steps involved in calculating the scattered electric field for a given 
orientation of the particle; then, a fifth and sixth step are required for averaging over orientation. The 
organization of these six steps into an EBCM algorithm is shown in Figure 3. We developed this 
logical structure specifically to take advantage of the ‘array processing’ or ‘vectorization’ features of 
modern computers. (Figure 3 refers to axisymmetric particles only; for non-axisymmetric particles, 
important sections of this logical structure would have to be re-formulated.) 

Step 1. Determination of n max in nose-on orientation 

The length N of the expansions of the incident and scattered fields in vector spherical harmonics (Eq. 
4) is not known a priori; it depends on the shape, size, and refractive index of the particle. In general, 
one should guess an initial N value, solve the scattering problem, and then increase N in steps until 
‘satisfactory’ convergence of the scattered field is obtained. This may, however, result in a prohibitive- 
ly lengthy procedure. 

Fortunately, for axisymmetric particles (see Sec. IV), it is possible to estimate the length of the expan- 
sions in a much faster way. For nose-on incidence (defined as 0 p = 0), by virtue of the rotational sym- 
metry of the particle about the direction of incidence, the only contribution to the scattered field 
comes from the single azimuthal mode m = 1 (as for a sphere). Therefore, it is advantageous to utilize 
nose-on orientation to determine the length of the n-expansion sufficient for convergence. Then we 
merely use the same upper limit n max for all the orientations. 

By so doing, we assume that n max is independent of the orientation of the particle — which is not true 
in general. For different orientations of the particle, deep downward spikes in the angular scattering 
pattern (representing a lot of cancellation among the terms in the spherical harmonic expansion) often 
require higher n max values in order to converge. However, we usually ignored this problem; our in- 
terest was mainly in randomly-oriented particles, where the scattering pattern is rather smooth — the 
spikes wash out in the orientation-averaging process. (For further discussion of this spike- washing-out 
process, but for spheres, see Sec. VIII.) 

We are not going to describe here the numerical procedures in the nose-on case, since they follow 
closely steps 2-4 in Fig. 3, which are described below. The only difference from the general case is that 
now the loop over the rn max disappears (m max must be 1), and only the orientation (0 p = 0, 4> p = 0 ) is 
considered. 

By solving the scattering problem for successively larger values of n max in the nose-on case, and stop- 
ping when the changes become ‘small enough’, n max is determined. Section VI (convergence pro- 
cedures) will flesh out the details of this operation. 
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(STEP 1) NOSE-ON (0 p = 0) 



find n max 


LOOP 

m max = 0) > n max 


(STEP 2) 

l,J,K,L; r,J’,K’,L’ 

(vec. 0 S ) 

(STEP 3) 

T m sub-matrices 

(vec. n) 

LOOP 

0 p : from 1 to N e 


LOOP 

4> p : from 1 to N^ 


(STEP 4) 

0 b> 4>b 

(vec. 0) 


a '/t’ b '/i 

(vec. n) 



(vec. n) 


F b q ( 0 b^b) 

(vec. 0) 


m max 



E F q (m, ( 0 ) 

(vec. 0) 


m = 0 



iq(0,m ma x ) 

(vec. 0) 


ff sca 

(vec. n) 


ff ext 


END 

4> p LOOP 


END 

0p LOOP 


(STEP 5) 

<iq( 0 )> 

(vec. B p ,4> p ) 

(STEP 6) 

°sca ^ > 4 °ext ^ 

(vec. n) 

IF 

(convergence over rn max ) EXIT 


END 

m max L00P 



Figure 3. Logical structure of EBCM algorithm, (vec. x) means that the indicated step is vectorized 
over the variable x. Other symbols are defined in the text. 
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Because of the splitting into (n max + 1) sub-problems (see Sec. IV), steps 2-6 in Fig. 3 appear inside a 
loop over the azimuthal index ‘m\ It is advantageous to check convergence of the series (32) at each 
step- of this m-loop, since convergence usually occurs well before the theoretical upper limit n max . 
Hence partial sums from m = 0 up to the current value of m (called m max in Fig. 3) are calculated. 
When the newest terms being added to these sums are small compared to the values of the sums, this is 
called ‘convergence over m max ’, and the algorithm is finished. Details of this procedure may be found 
in Sec. VI. 

Step 2: Calculation of the surface integrals 

Rapid calculation of surface integrals (12-18) is important, since otherwise an inordinate amount of 
computer time will be wasted in assessing convergence over n max (Step 1), as we shall see in Sec. VI. 
Hence, we used a Gaussian quadrature over the variable cos 0 S , which requires about a factor of 10 
fewer integrand evaluations than, for example, Simpson’s or Bode’s Rule (the latter was used in the 
EBCM code we received from Prof. Barber). Of course, Gauss rules assume a smooth particle shape, 
with no corners or edges; cubes and such-like are excluded. 

We have found that, for the Chebyshev particles, the number of Gaussian quadrature points N G must 
not be substantially smaller than n max ; otherwise, the surface integrals are not calculated accurately 
enough, and this may even cause the scattered field expansion to diverge. On the other hand, we have 
noticed that very high N 0 values (N G > > n max ) may also sabotage convergence, presumably because 
the error term R in the Gauss rule 


+ 1 N g 

( f(x)dx= £ Wj f(xj) + R(N g ) (34) 

-1 i=l 

is proportional to a high-order derivative of the integrand (which is not true in Simpson’s or Bode’s 
Rule): 


R(N g ) = 


2 2n + 1 (n!) 4 
(2n + 1) [(2n)!p 


f (2n) (x 0 ) 


(35) 


(Xj and W; are the absicissae and the weights, respectively, of the Gauss rule, and Xq is a point in the in- 
terval of integration). Examples of inaccuracy in Gauss rales for overly-large numbers of quadrature 
points can be found in Davis and Rabinowitz (1975). 


For the above-mentioned reasons, we have chosen to take N G = n max throughout the calculations, 
even though this choice wastes computer time during the determination of n max (because of the need to 
re-calculate surface integrals at each step of the iteration instead of just using the old values). 
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Step 3: Evaluation of the transition matrix 

Once the surface integrals have been calculated, one is in a position to evaluate the transition sub- 
matrices T m (33) in the body frame. This requires the inversion of the block-diagonal (axisymmetric- 
particle) form of the second matrix in (10), which is accomplished using the Gauss-Jordan method. 
We used Barber’s routine for this, rather than a library subroutine from LINPACK, because it takes 
advantage of the large number of known zeroes in the matrix. 

Barber (1977) discusses the ill-conditioning in this matrix inversion when any of the variables which 
define the scattering object (its size parameter, deviation from sphericity, and refractive index) is in- 
creased too far. This problem is a symptom that the transition matrix T (and therefore the scattered 
field coefficient vector) has become very sensitive to the length of the spherical harmonic expansions, 
so that no convergence can be determined. 

In the present study, we fixed the Chebyshev deformation parameter e at various values between 
-0.20 and +0.20; then, for each e, we increased the size of the particle until ill-conditioning was en- 
countered and convergence could no longer be achieved. Except for the very smallest size parameters, 
it was impossible to achieve convergence for magnitudes of e significantly above 0.20. 

Step 4. Evaluation of the scattered field for any orientation 

We can now solve the scattering problem for any orientation (0 p , «£ p ) of the particle. First of all, using 
(24) one must calculate the two arrays of angular coordinates (0 b , <p b ) which define the various scatter- 
ing directions s in the body frame (see Fig. 2). 

Then, the expansion coefficients of the incident field, (a^ 1 and b^ 1 ) are evaluated in the body frame us- 
ing (20). They are multiplied by the transition sub-matrix T m (30) to obtain the expansion coefficients 
of the scattered field for the considered azimuthal mode m (a M s and b^) in the body frame. 

Next, the scattered field coefficients are used in (23) to obtain the components of the far-field 
amplitude in the body frame; and then, these components are transformed back to the laboratory 
frame (Eqs. 26). This is done for a whole array of scattering angles 0 at once, in vector loops, as in- 
dicated on the right-hand side of Fig. 3. 

The computations carried out so far (Steps 2-4) all refer to a specific azimuthal mode m in the m max 
loop in Fig. 3. One must also update the running sums of the components F q ( m ) of the far-field 
amplitude vector corresponding to the partial fields E m s : 

m max 

F,(e,m„„)=E F,W(©) 

m = 0 


( 36 ) 
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Finally, the intensities i q and extinction cross-sections cr ext are evalulated from these running sums us- 
ing (28) and (29b). For the computation of the scattering cross-section <r sca (29a), all the a/ and b„ s 
calculated so far must be considered; thus, in (29a), stakes into account all possible combinations of 
a, n and m for n < n max and m< m raax . 

When studying scattering just in fixed particle orientation, convergence over m max should be checked 
at the conclusion of Step 4. On the other hand, when our interest lies in randomly-oriented particles, 
this test must be postponed until after the average over orientation is made (Steps 5 and 6). 


Step 5. Avering over orientation: intensities 

Orientation-averaging is accomplished using a two-dimensional Gaussian quadrature over orientation 
angles (0 p , </> p ) of the following generic form: 


<Y> = ~T~ j 

47T J o 


2tt 


f +1 

d<j> J Y(0, <f>) d(cosQ) 


1 

4-n: 


n 9 


i= 1 


£ w i Wj Y(Oj. </>;) 
j = i 


(37) 


(This is actually just the product of two one-dimensional Gauss rules, not a general two-dimension 
rule, for which no exact solution for the weights and abscissae is known.) Y is any real scattering 
quantity, such as intensity or cross-section (averaging the complex scattering amplitude would 
preserve phase relations, which we assume does not occur). cos(0;) and w ; are the N e abscissae and 
weights of the Gauss rule over (—1,1); and </>j and Wj are the N^ abscissae and weights of the Gauss 
rale over (0,2 it). 


In our computations, we have used a number of zenith orientations which depends on the particle size 
and on its deformation from a sphere: 


N e = (0.5x + 3) (20 |e| + l)/3 (38) 

For each zenith orientation 0 p i , the number of azimuth orientations is taken proportional to N e and 
to the area of the spherical zone at that zenith angle: 


M 0 = 2 N e sin 0 p>i 


(39) 


Both (38) and (39) are purely empirical rules which we developed. There is no a priori reason to expect 
that they will be satisfactory for shapes very different from Chebyshev particles. 

With only 750,000 words of memory on the CRAY-1, we were forced to set an upper bound of 10 on 
N e (we also never let it fall below 3). This is due to the necessity of storing two complex matrices, of 


dimension (no. scattering angles) x (no. orientations), containing the components of the vector far- 
field amplitude. This limits the total number of orientations to about 150, which may not be enough 
for large sizes and/or deformations. However, for axisymmetric particles, some simplifications enable 
us to double or quadruple the number of orientations. 



l)0 p = 90° (E° x >> = 0, E° y » = 1) 


2)4> p = 180° (E° x ,. = 1 , E° y . = 0) 

This simplification was found by Barber (personal communication). 
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Once the azimuthally-averaged cross-sections have been obtained by this trick, the orientation - 
averaged cross-sections follow directly from a Gauss quadrature over zenith orientation, as usual. If 
the particle has a plane of symmetry perpendicular to its rotation axis, this quadrature can be 
simplified as described in Step 5. 

Calculation of associated Legendre and Bessel functions 

In steps 2-4, it is necessary to calculate the associated Legendre function P n m (cos 0) with argument 
either cos 9 b (for obtaining the far-field amplitude) or cos O s (for obtaining the surface integrals). 
This is done by using forward recursion for fixed m, starting from n = m + 1 : 

P m n + i (cos©) = — - — - [(2n + 1) cos 0 P n m (cos0) + (n + rn)P m n _j (cos©)] (42) 

n— m 4" 1 

The recursion is initialized by: 

m 

P m m (cos©) = sin m 0 n (2k- 1) (43a) 

k= 1 

P ™ + 1 (COS0) = (2m + 1) cos0 P m m (cos0) (43b) 

In step 2, it is necessary to calculate the spherical Bessel functions of the first kind, j n , both for real (k 
r s ) and complex (m k r s ) arguments. It is also necessary to calculate the spherical Bessel function of the 
second kind y n for real argument k r s only, since it enters into the expression of the spherical Hankel 
function. 

For y n , we use a forward recursion, which is stable for real arguments (Lentz, 1975). However, the use 
of forward recursion for j n is unstable. For complex arguments, we calculated j n using the continued 
fraction method of Lentz (1975), while for real arguments we have used a variation on the Miller 
backward recursion scheme proposed by Ross (1972). 

In the original code provided to us by Prof. Barber, Bessel functions were generated, both for real and 
complex arguments, using the Miller backward recursion. We replaced such procedures because they 
can occasionally overflow for large arguments, 

Vectorization 

‘Vectorization’ means re-arranging an algorithm in order to make the longest loops the innermost 
ones. These innermost loops must also be free of recursions and other structures that would prevent 
their being executed all at once. 

We have deeply modified the original code we received from Prof. Barber in order to take advantage 
of the vector features of a CRAY-1 computer. We have vectorized any section of the code which 
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timing tests revealed to be using more than 10% of the computer time. Super-speed vector routines 
(furnished by CRAY) were also used for performing some matrix operations. The code is now about 
ten times faster than the original one that we received. Even so, this research has required on the order 

of 10 hours of CRAY time. 

Fig. 3 shows what loops were vectorized in the various steps of the. numerical computation. The nota- 
tion ‘vec. x 5 means that loops over the variable ‘x’ were made innermost, and vectorized. (If only ran- 
dom-orientation cross-sections are desired, Step 4 is replaced by a simplified version in which only 

loops over n are vectorized.) 

In the Bessel function routines, recursions prevent most of the loops over n from vectorizing. In the 
computation of the associated Legendre polynomials, either loops over 0 S or over 0 were vectorized, 
depending on the context (surface integrals or spherical harmonics expansion, respectively). 

Numerical cheeks 

Since considerable modifications of the code were made, we repeatedly checked the modified code 
against some results obtained with the original one supplied to us by Barber. We also checked the 
original code against some exact calculations of Asano and Yamamoto (1975) for spheroids in fixed 
orientation, for which agreement within 0.01% had been found (Mugnai and Wiscombe, 1980). 

After the modification of the code had been completed, we also checked it against exact calculations 
of Asano (1980) for randomly oriented spheroids. Very good agreement was still obtained (0.1%). 
However, we were not able to get convergence for elongations greater than 3:1 when the size 
parameter exceeded about 10. This is in line with a comment made by Waterman (1965) in his original 

paper on the EBCM: 

“One present drawback of the method is the poor numerical convergence of the truncation pro- 
cedure in dealing with more elongated shapes. This is of course not surprising, since one is de- 
parting from the nearly spherical shapes most ideally suited to the vector wave functions 
employed. As a rough rule-of-thumb, it appears that obstacles having a length-to-width ratio 
<2 can be treated quite well by truncation to at most 20-25 terms. ...” 

Holt (private communication) has informed us that other investigators who use the EBCM have been 
able to obtain convergence for more elongated spheroids. The explanation may lie in our stringent 
convergence procedures (see next section), compared to the relatively lax definitions of convergence 
used by others. However, we did not concern ourselves with this difficulty, because the many sensitivi- 
ty studies we carried out for Chebyshev particles convinced us that our numerical procedures are quite 
suitable for these particles. 
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In the EBCM, infinite series are truncated after a finite number of terms; integrals are numerically 
evaluated; and finite series may be truncated before their theoretical upper limit. The resulting solu- 
tion may be more or less ‘converged’, depending on our skill (there is more than a little art involved as 
well) in selecting the following three parameters controlling convergence: 

® the highest n- value, n max , used in the spherical harmonic expansions; 

® the number of Gauss quadrature points, N G , for performing the surface integrals; and 

• the number of different orientations of the particle, N or , used for calculating an orientation 
average. 

In addition, it is advantageous to check convergence on a fourth parameter: the highest m-value, 
rn max , used in the spherical harmonic expansions (this is not strictly necessary, since we know that, 
theoretically, m max = n max ). 

Our procedures for determining convergence over N G and n max have evolved from those used by 
Barber (private communication). We therefore begin by summarizing his method. 

Barber looks only at nose-on orientation of the particle. He guesses initial values for n max and N G . N G 
must be about 3-5 times n max , depending on the deformation of the particle from a sphere (M G must be 
large because Barber uses Bode’s Rule for calculating surface integrals). Then, n max is increased in unit 
steps, holding N G fixed, until ‘convergence’ of the scattered intensities is reached. 

‘Convergence’ is defined in the Cauchy sense: it occurs when two successive iterates are ‘close’. For 
Barber, ‘close’ meant that the percentage difference between the scattered intensities for n max and for 
( n max “ 1) was l ess than 1% at every scattering angle. (Clearly, the exact convergence he reached 
depended on the set of scattering angles he chose.) 

If convergence is slow, the Cauchy criterion can of course be misleading. Two successive iterates can 
be close, yet both can be far from the final solution. In practice, such pathological situations manifest 
themselves as a painfully slow downward trend in the inter-iterate difference. 

If convergence didn’t occur. Barber fixed n max at the highest value it had attained, and iterated over 
N g instead (by trial and error). If convergence over N G was attained, iteration over n max was re-started 
from its original (low) value, holding N G constant at its converged value. 
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We have modified Barber’s procedure for the following reason. Barber did EBCM calculations on a 
case-by-case basis, checking and tuning the convergence by hand until the results were satisfactory. 
Our goal was to study a great many cases, precluding the possibility of examining each one individual- 
ly; therefore, we had to invent an automatic procedure. Of course, we could simply have enforced a 
very stringent Cauchy criterion (0.01%, say); but this would have allowed convergence in only the 
easiest cases (those with the smallest deviation from sphericity). Furthermore, since our interest was in 
random orientation, we did not care to enforce a stringent convergence at every particular orientation 
and for every scattering angle. 


Convergence over M G 

In order to study convergence over N G , we carried out many tests for Chebyshev particles in nose-on 
orientation. In these tests, M G was either held fixed, or increased in direct proportion to n max : 

N G = k 0 n max. C 44 ) 


k 0 was varied between 1 and about 3. 

We found that, independent of n max , convergence over N c was reached for 1 < k 0 < 1.5 (larger k 0 
values were necessary when n max was below its convergence value). In general, we found that N G must 
be neither smaller, nor much bigger, than n max ; both extremes could sabotage convergence. 

Based on our tests, we decided to take k 0 = 1 throughout our calculations, for the following reasons: 

• increasing k 0 from 1 to 1.5 changes the intensities less than 1/10 as much as increasing the 
convergence value of n max by one; 

• this choice prevents the scattered field expansion from diverging; and 

• taking k 0 > 1 sometimes lowers the n max value for convergence, giving results of less ac- 
curacy than those obtained with k 0 = 1 . 

Eq. (44) is unfortunate in one sense: if N G were held constant, then each of the 8 matrices of surface 
integrals (e.g. Eq. 12) could be preserved for the next step of the iteration over n max ; only two new 
rows, and two new columns, would need to be added. Instead, each matrix has to be recalculated in 
toto. But at least (44) eliminates Barber’s trial-and-error double iteration (over N G and n max ) in favor 
of a single automatic iteration (over n max ). 


Convergence over n max 

The selection of n max Is critical in obtaining correct results, especially for fixed orientation of the parti- 
cle. Extreme caution is called for, first because bad results obtained with incorrect n max values are 
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difficult to detect; and second, because there are usually no independent methods to check the results 
(unless the particles are cylinders or spheroids). 

We adopted Barber’s suggestion to increase n max in unit steps during convergence iteration. However, 
we found that his convergence criterion (1% at all the angles) failed when the angular scattering pat- 
tern developed deep down-spikes, signalling cancellation among the terms of the spherical harmonics 
expansion at the spike angles. Convergence at those angles could not usually be achieved at all. 

However, our main interest was in random orientation, in which case the deep down-spikes one sees 
in fixed orientation are washed out. Therefore, we developed a convergence criterion with the follow- 
ing requirements in mind: 


• convergence should not be tested in the spikes 

® convergence should be summarized in a single number, representing some sort of average 
over all angles. 


We met the first requirement by considering, not all the N sca scattering angles, but only 90% of them, 
i.e. only 


Nsca = 0.9 N sca (45) 

scattering angles. We throw out that 10% of the angles where the percentage differences from one 
iteration to the next are largest. This is always sufficient to get rid of the spike situations. 

The second requirement is met by defining two numbers, <5 q (q = l or r), which are the root mean 
square relative differences between the intensities when n max is increased by one (but only at the 90% 
of the angles defined above): 


N, 


i = 1 


(®i> ^max) *q (®i> ^max 

i q ( 0 i> n max -l) 


2 j 'A 


‘Good’ convergence is then defined as 

S JL ^ ^min an d 5 r < S min 


(46) 


(47) 


where 5 min = 0.001 for our research. This corresponds to percentage differences in the intensities not 
exceeding about 0.1%, on the average. Thus we have relaxed Barber’s criterion in one sense — by 
disregarding 10% of the angles — and made it stricter in another sense — by demanding 0.1% rather 
than 1% convergence at the remaining angles. 

We have made many plots of the behavior of 5 q as a function of n max . Figure 4.1 shows a small but 
representative sample of this information; 5^ (dotted line) and <5 r (dashed line) are plotted vs. n max for 
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Figure 4,2. Parallel and perpendicular intensities i^and i r vs. scattering angle 0, for T 6 (0.10) 

particles in nose-on orientation with size parameter x = 10 and n max = 18, 30, 42, anc 
66 . (a-b) ijjfor 0-60 and 60-180 degrees respectively; (c-d) \ r for 0-60 and 60-180 degrees 
respectively. The solid line in each plot represents the “converged” result, 
corresponding to n max = 66. The various dashed lines correspond to lower values of 
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Figure 4.2 (Continued) 
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Figure 4.3. Parallel and perpendicular intensities i^and i r vs. scattering angle 0, for T 8 (0.10) 

particles in nose-on orientation with size parameter x = 8 and n max = 28,44, and 76: 
(a-b) i^for 0-60 and 60-180 degrees, respectively; (c-d) i r for 0-60 and 60-180 degrees, 
respectively. The solid line in each plot corresponds to the global minimum n max = 

44 in Figure 4.1(d). The dotted line corresponds to n max = 28, the dashed line to n max 
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Figure 4.3 (Continued) 





Figure 4.4. Parallel intensity i^vs. scattering angle 0, for T 6 (0.10) particles in random orientation 
with size parameter x = 10 and n max = 18,30,42“ and 86: (a) 0-60 degrees; (b) 80-180 
degrees. The solid line and the various dashed lines have the same meaning as in 

Figure 4.2 
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four different Chebyshev particles. Fig. 4.1(a) provides an example of how convergence becomes more 
difficult when the size of the particle increases; Fig. 4.1(b), of how it becomes much more difficult 
when the deformation from sphericity increases. 

Figure 4.1 shows clearly that <5 q does not decrease monotically when n max is increased, but rather has 
oscillations superimposed on its average downward trend. These oscillations have a period equal to the 
order n of the Chebyshev polynomial T n defining the particle! This surprising (and so far unexplained) 
result explains why convergence over n max is very slow for high Chevyshev orders. 

Figure 4.1 also shows that 5 q eventually stops decreasing, and may actually start increasing again, as 
n max continues to increase. This behavior may be due to 

• accumulation of round-off errors; or 

• the possibility that the spherical harmonics expansion, at least for concave particles, is an 
asymptotic rather than an exact series (this may be a consequence of the Rayleigh hypothesis 
which we have made); it is well known that there is an optimal number of terms to keep in an 
asymptotic series, beyond which it becomes less and less accurate. 


These plots show explicitly how our definition of ‘convergence’ may fail: 5 q may ‘bottom out’ before 
(47) is satisfied. However, our convergence studies also revealed that satisfactory convergence of the 
scattered intensities can often be achieved even for n max values violating (47). 


For this reason, we have introduced a second, less stringent convergence criterion based on the 
quantity 



_ Zi 


Vi (5j + 62) 


(48) 


We stop the n max iteration before ‘good’ convergence is reached, if 6^ r reaches a minimum lower than 
10 times 6 mjn ( = 0.01 in our study), then starts increasing again . We call this convergence ‘fair’. It 
produces accurate enough results for a randomly oriented particle. However, it may give unacceptable 
errors for fixed orientation. 

By introducing ‘fair’ convergence, it was possible to obtain results for many more cases than ‘good’ 
convergence would have allowed. It also shortens the iteration over n max (which may otherwise be in- 
ordinately lengthy due to the oscillations in the 5 q vs. n max curves). This in turn reduces computer time 
spent in orientation-averaging, which goes roughly as n max to the 2.5 power. 

We have singled out the two particles in Fig. 4.1(c-d) for particular consideration, since they are 
representative of cases of slow or difficult convergence, respectively. The first is a T 6 particle with x = 
10; the second, a T g particle with x = 8. Both particles have deformation e = 0.10. 
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Parallel and perpendicular intensities for these two particles are shown in Figs. 4.2 and 4.3, respective- 
ly, for various n max values and for nose-on orientation. (Studies we made for oblique orientation did 
not differ qualitatively.) There are separate plots for 0-60 and 60-180 degrees in order to gain resolu- 
tion. 

For the T 6 ( + ) particle, Fig. 4.1(c) shows that a deep local minimum, both for 6^and o r , is reached 
already for n max = 18; ‘fair’ convergence is reached for n max = 30; and ‘good’ convergence, for n max 
= 42. For n max = 66, c>£ and 5 r are even lower yet. The solid curve in each plot in Fig. 4.2 represents 
the ‘converged’ result, corresponding to n max = 66. The various dashed curves correspond to lower 
values of n max . This enables us to watch the convergence process ‘in action’, so to speak. 

The trend towards convergence in Fig. 4.2 is quite evident. Moreover, we find that: 

® the main characteristics of the scattering pattern are already produced by n max = 18 

® perpendicular intensities converge much more rapidly than parallel ones 

• ‘fair’ convergence does not produce accurate enough results for fixed orientation 

® ‘good’ convergence produces almost perfectly converged intensities 

The T g particle is truly a marginal case; even a small increase in its size or deformation causes con- 
vergence to fail entirely. As it is, Fig. 4.1(d) shows that ‘fair’ convergence is reached for n max = 28, 
while ‘good’ convergence is not achieved at all. But there is a global minimum at n max = 44. Beyond 
that, d^and 5 r increase, but are still less than 0.01 in two local minima at n max = 60 and 76. Therefore, 
in Fig. 4.3 we have plotted parallel and perpendicular intensities for n max = 28, 44, and 76. We find 
that: 


• intensities for n max = 28 (dotted line) and 44 (solid line) are quite similar 

• intensities for n max = 76 (dashed line) are rather different from the first two. 

This T 8 example clearly shows that taking too many terms in the spherical harmonics expansion is just 
as deadly as taking too few. Therefore, it is unwise just to arbitrarily choose a high value of n max 
without checking convergence. For the same reason, even the first guess for n max in the convergence 
iteration shouldn’t be too high. But we have found that, even for ‘fair’ convergence, n max is always 
bigger than the number of terms, 


n sph = x + 4x' /3 + 2 (49) 

required for convergence of the scattered field from an equal-volume sphere (Wiscombe, 1980). 
Therefore, we suggest utilizing n max = n sph as a first guess. 
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Table 1 shows the n max and devalues corresponding to a selection of our Chebyshev particles, as well 
as what kind of convergence (‘fair’ or ‘good’) was obtained. 

Figures 4.4 and 4.5 are analogous to Figs. 4.2 and 4.3, respectively, but for random orientation. Only 
parallel intensities are shown, because the curves for various n max values are virtually indistinguishable 
for perpendicular intensity. Hence, for checking random-orientation convergence, it may be sufficient 
to look only at parallel intensity, although we have not done so in our calculations. 

Table 1. Convergence parameters for a selection of Chebyshev particles. F and G refer to ‘fair’ and 
‘good’ convergence, respectively. The other parameters are defined in the text. 


particle 

type 

n max 

^max *Vph 

Hr 

m max 

N or 

T 3 (0.05) 

G 

14 

1 

5e-4 

9 

48 

T 3 (0.10) 

G 

24 

11 

7e-4 

9 

132 

T 3 (0.15) 

F 

37 

24 

3e-3 

9 

132 

T 2 (±0.20) 

G 

17 

4 

5e-4 

8 

400 

T 4 (±0.05) 

G 

13 

0 

2e-4 

8 

56 

T 4 (±Q.15) 

F 

21 

8 

7e-3 

9 

284 

T 6 (±0.05) 

G 

13 

0 

3e-4 

8 

56 

Tg(±0.05) 

G 

16 

3 

9e-4 

8 

56 

T 8 (±0.10) 

F 

20 

7 

8e-3 

8 

144 

T 3 (0.05) 

F 

20 

0 

2e-3 

15 

132 

T 3 (0.!Q) 

G 

35 

15 

4e-4 

15 

132 

T 3 (0.15) 

F 

43 

23 

9e-3 

15 

132 

T 2 (-0.20) 

G 

28 

8 

4e-4 

15 

584 

T 2 ( + 0.20) 

F 

24 

4 

5e-3 

12 

584 

T 4 (±0.05) 

G 

20 

0 

3e-4 

14 

144 

T 4 (-0.15) 

F 

54 

34 

le-2 

13 

584 

T 6 (±0.05) 

F 

21 

1 

4e-3 

14 

144 

T 8 (±0.05) 

F 

22 

2 

3e-3 

14 

144 

T 8 (±0.10) 

F 

46 

26 

9e-3 

14 

372 

T 3 (0.05) 

F 

27 

0 

3e-3 

21 

132 

T 3 (0.10) 

F 

39 

12 

4e-3 

21 

132 

T 2 (±0.20) 

G 

41 

14 

9e-4 

20 

584 

T 4 ( ± 0.05) 

G 

27 

0 

6e-4 

20 

284 

T 6 (±0.05) 

F 

28 

1 

6e-3 

20 

284 

T 3 (0.05) 

F 

34 

1 

3e-3 

27 

132 

T 3 (0.10) 

F 

49 

16 

3e-3 

27 

132 

T 2 ( ± 0.20) 

F 

47 

14 

9e-3 

24 

584 

>ri 

© 

© 

+1 

G 

34 

1 

8e-4 

26 

372 
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There are several important observations to make about Figs. 4.4 and 4.5: 

• the spread between the curves for various n max values is an order of magnitude less than for 
nose-on orientation 

« ‘fair’ convergence is good enough for random orientation; the worst errors are below 2% or 

so 

• even the n max = 18 intensities in Fig. 4.4 are within 10% of the converged values 

9 the n max = 28 and 44 curves in Fig. 4.5 are virtually indistinguishable, whereas the n max = 76 

curve departs significantly from them. 

This last observation (which was also noted in nose-on orientation) again shows the danger of taking 
too many terms in the spherical harmonic expansion. 


Convergence over m max 

We have found that, for a given n max , convergence over m is usually reached for m max ^ n max - ^ i® 
therefore advantageous to introduce a convergence criterion over m max in order to shorten the com- 
putation. 


We adopted, an m max convergence criterion similar to ‘good’ convergence for n max (46-47), except that 
all the scattering angles are now taken into account, because down-spikes do not seem to occur in the 
m-sum (apparently there is never the kind of massive cancellation of terms that occurs in the n-sum). 
For a particle in fixed orientation, the analogue to (46) is therefore: 


[ ^sca 

j i y 

iq (9i> ^max) iq (®i> *Wnax 

2 

N sca ", 

iq ^max ~^) 



(50) 


where 


iq (®> 'Wnax) 


m max 

£ F q (m) (0) 

m = 0 


(51) 


and F q ( m ) is the q-component of the far-field amplitude of the partial field E m s . 

For a randomly-oriented particle, convergence over m max is checked by making use of an exactly 
analogous criterion, except that i q is replaced by <i q > (see 37). 

% (dotted line) and <5 r (dashed line) are plotted vs. m max in Figure 4.6, for two different Chebyshev par- 
ticles in random orientation. There is a striking trend towards convergence as soon as these quantities 
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become smaller than about 0.001; hence, convergence over m max was assumed to have taken place as 
soon as 8% and <5 r fell below this limit. 

When only random-orientation cross-sections are desired, and the short-cut described in ‘Step 6" of 
Sec. V is employed, m max convergence is based on the following quantity: 


1 

N e 

v 1 

a <t> (®p, i’ m max) °<j> (®p, i> 

2 

N e 

Ld 
i = 1 

®<t> (®p, i» m max 



(52) 


where <j> means one of the two azimuthal orientations (4> p = 90, 180 degrees) used for computing the 
random-orientation cross-section (Sec. V), and N e is the number of zenith orientations. In general, 
this leads to the same m max values as applying the 5 q criterion. 

Table 1 shows m max values for various particles. They range between about 14 and 3 A of n max , showing 
that checking convergence over a m max reduces computer time by roughly a factor of two. 

Convergence over N e 

When we compute the scattered field from a randomly-oriented particle, the number of azimuth 
orientations N^ for any given zenith orientation is automatically assigned (39). Therefore, convergence 
is checked only over N e . 

Computer time and storage constraints have prevented us from testing convergence over orientation 
for many of the cases in this compendium. Ideally, computations should be done for successively 
larger N e values until the final results converge to a specified accuracy. But in that case, in order to 
keep computer time within reasonable bounds, one would need to store all the T m transition matrices; 
and this was beyond the capacity of the NCAR CRAY-1 (750,000 words). 

We therefore resorted to testing the random-orientation results in selected cases in which N 0 did not 
seem large enough. We temporarily modified the EBCM code to allow twice the usual maximum 
number of orientations (584 for T n with even n, 1 32 for odd n). For a T 2 ( - 0.20) particle with x = 20, 
and a T 3 (0.15) particle with x = 10, the largest changes due to doubling the number of orientations 
were about 1-2%. 

Suggested improvements in the convergence procedures 

We wish now to suggest how the convergence criteria could be improved. However, in spite of our ex- 
tensive experience, we are still unable to give general and completely automatic convergence pro- 
cedures. Our suggestions are probably most appropriate to the Chebyshev particles; they may however 
serve as guidelines for other shapes. 
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We have not introduced the following changes in our computer code yet. However, we believe that 
they could speed up the convergence procedures and at the same time produce even more accurate 

results. 

First, even though ‘fair’ convergence seems to produce accurate enough results for randomly-oriented 
Chebyshev particles, it might not be restrictive enough for other shapes. Besides, ‘fair’ convergence 
should always be avoided for fixed orientations. We suggest removal of the ‘fair’ convergence option 
and the raising of the value of 5 min instead. 5 min = 0.003 should produce accurate enough results for 
random orientation; and <5 min = 0.001 for fixed orientation. 

Second, we propose an improvement in the iteration over n max whenever the code discerns a periodici- 
ty in the variation of 8 % and 5 r with n max . The iteration over n max should then proceed in steps equal to 
the periodicity, rather than in steps of one. This method could bring a considerable speed enhance- 
ment for those particles with long periodicities (See Fig. 4.1c). If ‘good’ convergence does not occur in 
this case, one should use the results from the overall minimum of 5fl r , comparing them with results for 
the previous and following minima to determine if the convergence is reliable (see Figs. 4.1(d), 4.3 and 
4.5). 


Third, the convergence over number of orientations could be handled in a more economic fashion, by 
using ‘adaptive’ Gaussian quadrature (Patterson, 1973). Here, one uses values of N e in the sequence 
1, 3, 7, 15, 31, 63, and so on. Each successive Gauss rule employs all the points used by its predecessor. 
Hence, one can select a high value (say 15) of N e , and obtain orientation averages for the lower values 
(say 3, 7) for free. Comparing the various values will show how the convergence is going. 

Finally, we propose yet another way to test convergence over number of orientations. We have noticed 
that backscattering at 180 degrees is very sensitive to number of orientations. Therefore, one could run 
the EBCM just for a scattering angle of 180 degrees, testing convergence for successively increasing 
values of N e ( perhaps as described in the previous paragraph). This would be relatively cheap, com- 
putationally, since computer time for calculating random-orientation intensities is roughly propor- 
tional to the number of scattering angles. Once convergence was achieved, the EBCM could be re-run 
with the full complement of scattering angles, holding N e fixed at its converged value. In this pro- 
cedure, it would be desirable to save all the T m transition matrices to avoid recalculation. 



VII. DEFINITIONS OF QUANTITIES TO BE PLOTTED 


In order to deal with numbers of order unity, it is customary to define the efficiency factors 

Qsca = °sca /7r r ev 2 (53a) 

Qext — ^ext^ 71 " r ev 2 (53b) 

in which the cross-sections are normalized by the projected area of the equal-volume sphere. The 
single-scattering albedo is then defined as 

w — Qsca^Qext (54) 

The unpolarized phase function (which is proportional to the ‘differential scattering cross-section’ 
used by the inventors of the EBCM) is defined as 

P(cos 9) = 2 (ijj + i r )/(x 2 Q sca ) ( 55 ) 

where i^and i r are the parallel and perpendicular scattered intensities defined in (28). P, being essen- 
tially a probability of scattering at angle 9 from the direction of the incident unpolarized radiation, is 
normalized to unity: 


(Vi ) j P(cos 9) d(cos 9) = 1 

-i 


( 56 ) 


In practice, we only calculate intensities at discrete angles; therefore, this normalization relation can 
test whether or not we have picked enough angles. With our usual 1 1 1 angles, and a Simpson Rule for 
the left-hand side of (56), we always obtained numbers in the range 0.997 to 1 .003, indicating that our 
Simpson-Rule calculation of the asymmetry factor 


g = 



cos 9 P(cos 9) d(cos 9) 


(57) 


is probably accurate to 2 decimal places, with an angular-quadrature error of perhaps ± 1 in the 3rd 
place. 

We use a trick to calculate unpolarized phase functions in the EBCM. The straightforward way would 
be to calculate ig and i r first for incident parallel-polarized light, then for incident perpendicular- 
polarized light, and add them up. (For non-spheres, even in the random orientation, there is always 
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cross-polarization, meaning that incident parallel or perpendicular light produces non-zero values of 
both L and i r ) Instead, we assume incident light with linear polarization at 45 degrees to the scattering 
plane. This gives the same answer for i^and i r (but not for the remaining two Stokes parameters), with 
only half the computation. This trick works because, in our cases, the 4x4 Mueller matrix has 
vanishing 2x2 sub-matrices in the upper right and lower left corners: nose-on orientation falls into case 
2 of van de Hulst (1981), sec. 5.22; random orientation falls into case 6 of that same section. 


Clearly, for a general non-spherical particle in fixed orientation, the phase function would depend on 
the azimuth angle between the scattering plane (variable) and the plane containing both the incident 
direction and the reference axis (fixed), as well as upon the scattering angle 9. However, for the par- 
ticular situations we consider — either nose-on (incident radiation along the particle’s axis of rotation) 
or random orientation — this azimuthal dependence disappears. To show the phase function for ar- 
bitrary particle orientation would really require a 3-D plot. 

The final quantity of interest is the degree of polarization: 

d(9) = (i r - fg)/ (i r + ig) (58) 

Like the efficiency factors, this is of order unity, and in fact is bounded between -1 and + 1. 

Experimenters prefer looking at phase function and degree of polarization because both are indepen- 
dent of the absolute calibration of the instruments (provided the calibration is the same for i^ and i r ). 
d(6) is particularly nice because, unlike the phase function, it can be determined accurately without 
knowing the scattering cross-section. By reading both phase function and degree of polarization from 
our plots, and looking up Q sca in Table 3, the reader can recover i^ and i r . This is burdensome, we 
know, but this compendium would have been too long if phase function, degree of polarization, and 
Intensities were all plotted. 

Scattering cross-section requires an integration over all scattered intensities. Unfortunately, in the 
laboratory, it is nearly impossible to position the detector at scattering angles 0 to 5 degrees, because 
then it receives incident as well as scattered radiation. Angles 175-180 degrees are generally 
unreachable as well, since here the detector either blocks the incident beam, or is blocked by the device 
emitting it. The loss of 0-5 degree information is particularly harmful, because this is in the diffraction 
region, where intensities are largest. 


VIII. SIZE-AVERAGING OF SPHERICAL RESULTS 


After examining many hundreds of plots comparing spherical to non-spherical scattering results, we 
came to feel that comparing to a sphere for a specific value of size parameter ‘x’, far from being the 
most honest comparison, was in fact the most misleading. Often the pattern of spherical-nonspherical 
differences would alter radically with only a slight change in x; and invariably that alteration was caus- 
ed by the sphere. The random-orientation nonspherical results were rock-steady by comparison. 

The problem lies, of course, in the multiple periodicities and spikes in the Mie results as a function of 
x. This is collectively known as the “ripple structure”. It is very special to a sphere, and is caused by 
the interference between surface waves and the internally transmitted and reflected radiation. Ripple 
tends to be damped out by nonsphericity, and especially by nonsphericity plus averaging over orienta- 
tion. (However, a high imaginary index will damp the ‘short-cutting’ surface waves, and hence the 
ripple, even for a sphere (van de Hulst, 1981).) 

The fastest oscillation in Mie quantities has a period of about 0.8 in x for refractive indices typical of 
aerosols in the shortwave spectrum (Dave, 1969)'. This has recently been explained as a ‘forward glory 
oscillation’ by Nussenzveig and Wiscombe (1980), who give an analytic formula for the period as a 
function of real refractive index. For m re = 1.5, this period is 0.7393 for spheres with size parameter 

x> > 1. 

To see the extent to which this period applied for smaller values of x, we made a number of high- 
resolution plots of Mie single-scattering albedo, w, and asymmetry factor, g, as a function of x, using 
the Mie algorithms of Wiscombe (1979, 1980). The complete course Of these quantities, from x = 0 to 
30, is shown in Fig. 5(a); then, in Figs. 5(b-i), the sub-ranges x = 2-6, 6-10, 10-15 and 15-20 are shown 
in detail (solid lines). For w below x = 10, the peak-to-peak period is about 0.85, while the trough-to- 
trough period is in the range 0.71 to 0.76. Both periods settle down to around 0.74 ±0.02 for x > 10 . 
For g, the period is consistently around 0.78 for all x <25. 

Since the deviations from the theoretical period are relatively small, we averaged the Mie single-scat- 
tering albedo and asymmetry factor over an interval 5x = 0.7393, using 201 equally-spaced points 
centered on each value of x. These averages are tabulated for integer values of x in Table 2 alongside the 
exact (unaveraged) values. Then, they are plotted as dotted lines in Figs. 5(b-i). The resulting curves are 
almost completely smooth, with only a faint residual oscillation due to not averaging over the exact 
period. 

The averaging was actually done over cross-sections and intensities (which makes more physical sense) 
rather than directly on o> and g: 


x2 Qsca /x2 
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MIE SIZE PARAMETER 


Figure 5. Mie single scattering albedo « and Mie asymmetry factor g vs. Mie size parameter x: 

(a) for the full range x = 0-30 considered in this study; (b-e) exact (solid line) and size- 
averaged (dotted line) oj in the sub-ranges x = 2-6, 6-10, 10-15, and 15-20; (f-i) exact (solid 
line) and size-averaged (dotted line) g in the same sub-ranges. 
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Qext = x2 Qext/X 2 

(59b) 

w avg — Qsca^Qext 

(59c) 

§avg = X 2 gQsca/x 2 Q SC a 

(59d) 


The overbar indicates size-averaging over a rectangular size distribution of spheres spanning the range 
Ax. 

Turning now to the phase function, we observe that the deep downward spikes in exact Mie phase 
functions, separating the (approximately) ‘x’ peaks between 0 and 180 degrees, must be due to very 
particular cancellations of terms in the Mie series. Small changes in x cause big changes in the location 
of these spikes. Hence even a modest amount of size-averaging reduces them considerably. The peaks 
between the spikes, on the other hand, are less affected by size-averaging. Thus, the general effect of 
size-averaging is to smooth the Mie curves considerably. 

Figure 6 illustrates this phenomenon. The thin solid curve is a size-averaged (Ax = 1) Mie phase func- 
tion for x = 15 in the angular range 80-180 degrees; the thick solid curve, the exact Mie result. The 
dotted and dashed lines show the exact and size-averaged < T 4 (-0.10) > phase function; these bear 
out our remark above, that nonspherical results tend to be rock-steady under size-averaging, com- 
pared to spherical results. 

Away from scattering angle 0 = 0, Mie phase function values tend to have several incommensurate 
periods as a function of x — not just 0.8. For example, at 180 degrees (the worst angle), Shipley and 
Weinman (1979) find periodicities of 0.42, 0.81, 1.1, etc. using power-spectrum analysis. Thus, there 
is no clear theoretical justification for averaging over an interval Ax = 0.7393. In fact, when we used 
this fixed averaging interval, it gave too much averaging for small x, and not enough for large x. 
Hence, we found it more convenient to use Ax = O.lx for all size-averaged spherical phase functions. 
This of course aliases all the oscillations, but still leads to an acceptable degree of smoothing of the 
most unrepresentative parts of the Mie curves. 

To maintain consistency between our treatment of single-scattering albedo, asymmetry factor, and 
phase function, we will also compare (in Sec. IX) nonspherical w and g with their spherical counter- 
parts averaged over Ax = O.lx. That is why the latter values were also included in Table 2. Except for 
x < 3, this causes very little change in the spherical-nonspherical percentage difference plots. 


Goedecke, in a comment at the end of Wiscombe and Mugnai (1980), suggested using the pro- 
jected-area distribution of the non-spherical particles as the size distribution for averaging the 
spherical results. Since many of our particles exhibit 10% deviations from a sphere, we had something 
like this in mind in picking Ax = O.lx. However, in detail, Goedecke’s idea would seem more ap- 
plicable to large particles (x > > 1) than to the ones we studied, because projected area and scattering 
cross-section are closely related only in the geometric-optics limit. 
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SCATTERING ANGLE 


Figure 6. Exact (thick solid line) and size-averaged (thin solid line) Mie phase functions, and 

exact (dotted line) and size-averaged (dashed line) <T 4 (-0.10)> phase functions, for 
size parameter x = 15. Shows the relatively small effect of size-averaging on 
nonspherical results. 
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Since we will not be showing size-averaged scattered intensities, it is worth noting that the effect of 
size-averaging is much greater on i r than on i^. For i r , it is typical for down-spikes to be raised by fac- 
tors of 5 to 10, and sometimes even by factors of 100. For both intensities, the smoothing is usually 
greatest in the angular region 80-150 degrees. The oscillatory structure in the 0-80 and 150-180 degree 
regions is more resistant to size-averaging. 

Looking ahead to Figures 9 and 10, the reader can see examples of exact and Ax = O.lx size-averaged 
spherical phase functions and degrees of polarization. 


Table 2. Mie single-scattering albedo (gj) and asymmetry factor (g): exact, Ax = 0.7393 size-averaged 
and Ax = O.lx size-averaged, x is Mie size parameter. 


X 

to 

“avg 

(Ax = 0.74) 

“avg 

(Ax = O.lx) 

g 

Savg 

(Ax = 0.74) 

Savg 

(Ax = O.lx) 

1 

0.787 

0.824 

0.788 

0.200 

0.304 

0.203 

2 

0.907 

0.908 

0.907 

0.634 

0.640 

0.634 

3 

0.922 

0.919 

0.921 

0.747 

0.748 

0.747 

4 

0.910 

0.909 

0.909 

0.772 

0.772 

0.772 

5 

0.873 

0.877 

0.876 

0.751 

0.751 

0.752 

6 

0.819 

0.816 

0.816 

0.701 

0.686 

0.690 

7 

0.736 

0.729 

0.730 

0.611 

0.604 

0.604 

8 

0.678 

0.688 

0.689 

0.612 

0.631 

0.633 

9 

0.707 

0.712 

0.713 

0.748 

0.748 

0.749 

10 

0.742 

0.738 

0.737 

0.829 

0.827 

0.826 

12 

0.709 

0.708 

0.708 

0.860 

0.859 

0.859 

14 

0.638 

0.637 

0.639 

0.838 

0.837 

0.839 

15 

0.647 

0.645 

0.645 

0.847 

0.850 

0.852 

16 

0.660 

0.663 

0.663 

0.868 

0.871 

0.871 

18 

0.656 

0.654 

0.653 

0.884 

0.883 

0.883 

20 

0.600 

0.600 

0.603 

0.875 

0.876 

0.878 

22 

0.609 

0.608 

0.608 

0.900 

0.899 

0.899 

24 

0.612 

0.612 

0.609 

0.913 

0.913 

0.912 







IX. SINGLE-SCATTERING ALBEDO AND ASYMMETRY FACTOR 


Single-scattering albedo and asymmetry factor, aside from their intrinsic interest, play a crucial role in 
all simple models of multiple scattering (e.g., Joseph, et. al., 1976). Figure 7 shows percentage dif- 
ferences between spherical and nonspherical single-scattering albedo. Figure 8 is similar, but for the 
axymmetry factor. Results are shown only for the MIXTURE, and for T n particles up through T 4 ; we 
had much less data for T 6 , T 8 and T 20 particles, and it showed deviations of less than 1-2% from the 
spherical results. Because of computer time limitations, we could not collect enough data to make 
these curves smooth; however, we believe that they are not seriously aliased, even if they present a 
rather unattractive appearance. 

Percentage differences for the MIXTURE (Figs. 7a and 8a) are with respect to exact (dotted line), Ax 
= O.lx size-averaged (dashed line) and Ax = 0.7393 size-averaged (dot-dash line) spherical results. It 
doesn’t matter much which Mie curve is used, as the reader can easily see, and the differences from the 
sphere are generally below 3%. 

In Figs. 7(b-d) and 8(b-d), percentage differences for T n particles are only with respect to the Ax = 
O.lx size-averaged Mie results; on each of these plots, the dotted, dashed, and dot-dash lines represent 
the various | e | values shown, in increasing order. In the T 2 and T 4 cases, there are two of each kind of 
curve; one is for the positive value of e, the other for the negative value. The two curves are not 
separately identified because they are so close to one another. 

The spherical-nonspherical differences in Figs. 7 and 8 were usually much greater than the spread be- 
tween the exact and size-averaged spherical curves in Fig. 5. Also, spherical and nonspherical asym- 
metry factors almost always differed in the second decimal place, so the angular-quadrature error in 
computing nonspherical g’s (in the 3rd decimal place) is not significant either. This gives our conclu- 
sions about the effect of nonsphericity on w and g a certain robustness. 

The ‘noise level’ on the percentage differences in Figs. 7 and 8 is somewhere between 0.1 and 1%. The 
lower limit has to be 0. 1 % because we always rounded the nonspherical quantities to 3 significant 
digits. The upper limit is due to obtaining only ‘fair’ rather than ‘good’ EBCM convergence in some 
cases. 

Because we only present plots of spherical-nonspherical differences in « and g, and because plots of w 
and g are difficult to read to 3 significant digits, we have tabulated them, as well as Q sca and Q abs , for 
each Chebyshev particle studied, in Table 3. (Q sca is necessary to convert phase functions back to in- 
tensities; see Eq. 55.) The particular value of ‘x’ at which each tabular column ends generally 
represents the largest integer for which we were able to obtain at least ‘fair’ convergence in the EBCM. 
For e = 0.05 (and sometimes even for e = 0.10), however, calculations are stopped before this barrier 
is reached. 
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MIXTURE 



Figure 7. Percentage difference between spherical and nonspherical single-scattering albedo, 
vs. size parameter x: (a) MIXTURE (a uniform blend of all particle shapes considered); 
(b) T 3 for e = 0.05,0.10, and 0.15; (c) T 2 for e = -0.10, +0.10, and -0.20, +0.20;(d)T 4 
for e = 0.05, +0.05, -0.10, + 0.10, and -0.15, + 0.15. In (a), differences are with respect 
to exact (dotted line), Ax = O.lx size-averaged (dashed line), and Ax = 0.7393 size- 
averaged (dot-dash line) spherical results. In (b-d), differences are only with respect to 
the Ax =0.1x size-averaged spherical results; the dotted, dashed, and dot-dash lines 
represent the various e values shown, in increasing order. Dog-legs are due to large 
steps between computations; there may be finer structure not revealed by these plots. 
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Figure 7 (Continued) 
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MIXTURE 



(T 3 ~) l£l = 0. 05,0. 10,0. 15 



Figure 8. Percentage difference between spherical and nonspherical asymmetry factor, vs. size 
parameter. Same cases, and same meanings for the lines, as in Fig. 7. 
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<V),(V> |£|=0. 10,0.20 



<T 4 ->,(T 4 + > IeI=0.05, 0.10, 0.15 



Figure 8 (Continued) 




SINGLE SCATTERING FROM NONSPHERICAL CHEBYSHEV PARTICLES 


The set of data in Table 3 is one of the key products of this study. It is much more extensive than any 
which has been computed before; and the entire experimental literature does not contain even 1/10 
this much data on nonspherical co and g. Furthermore, together with Asano’s (1980) study for 
spheroids, this is the only study to show the effect of nonsphericity plus random orientation on both w 

and g simultaneously. 


It may be of some interest for the reader to compare these results with the predictions of the Pollack/ 
Cuzzi (1980) semi-empirical theory, which are that: 


^nonspher ' > ^spher 

(60a) 

Snonspher ^ Sspher 

(60b) 


where ‘spher’ refers to an equal-volume sphere, as in our case. 
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Table 3. Scattering efficiency Q sca , absorption efficiency Q abs , single-scattering albedo w, and 
asymmetry factor g for spheres (S), Ax = 0. lx size-averaged spheres (S avg ), and various 
Chebyshev particles as a function of size parameter x (equal-volume size parameter 
for the nonspherical particles). 


SCATTERING EFFICIENCY Q sca 


X 

S 

c 

°avg 

T 3 (0.05) 

T 3 (0.10) 

T 3 (0.15) 

T 2 (-0.10) 

T 2 (0.10) 

2 

1.66 

1.66 

1.65 

1.65 

1.63 

1.64 

1.64 

3 

3.05 

3.05 

3.05 

3.05 

3.04 

3.04 

3.04 

4 

3.54 

3.57 

3.54 

3.53 

3.53 

3.57 

3.59 

5 

3.25 

3.21 

3.24 

3.23 

3.22 

3.24 

3.26 

6 

2.34 

2.29 

2.35 

2.38 

2.42 

2.37 

2.36 

7 

1.51 

1.55 

1.52 

1.59 

1.73 

1.62 

1.62 

8 

1.35 

1.39 

1.38 

1.47 

1.61 

1.42 

1.44 

9 

1.65 

1.65 

— 

— 

— 

1.64 

1 .64 

10 

2.00 

1.96 

1.98 

1.99 

2.03 

1.92 

1.90 

12 

1.78 

1.79 

1.78 

1.75 

1.71 

1.79 

1.82 

14 

1.35 

1.36 

1.35 

1.36 


1.42 

1.42 

15 

1.41 

1.43 

1.43 

1.46 


1.45 

1.43 

16 

1.57 

1.57 

1.56 

1.59 


1.54 

1.51 

18 

1.57 

1.56 

1.56 

1.53 


1.53 

1.56 

20 

1.30 

1.31 

1.29 

1.30 


1.37 

1.38 


X 

T 2 (-0.2Q) 

T 2 (0.20) 

T 4 (-0.05) 

T 4 (0.05) 

T 4 (-0.1Q) 

T 4 (0.10) 

2 

1.61 

1.59 

1.65 

1.65 

1.64 

1 .64 

3 

3.00 

3.00 

3.05 

3.05 

3.04 

3.04 

4 

3.59 

3.65 

3.55 

3.55 

3.57 

3.57 

5 

3.26 

3.32 

3.25 

3.24 

3.26 

3.23 

6 

2.48 

2.54 

2.35 

2.35 

2.36 

2.38 

7 

1.82 

1.87 

1.54 

1.55 

1.60 

1.65 

8 

1.57 

1.55 

1.38 

1.37 

1.47 

1.45 

9 

— 

— 

— 

— 

1.76 

1.75 

10 

1.83 

1.75 

2.00 

2.01 

2.04 

2.06 

12 

1.76 

1.89 

1.80 

1.79 

1.86 

1.85 

14 

1.60 

1.62 

1.34 

1.35 

1.36 

1.36 

15 

1.58 

1.48 

1.39 

1.39 

1.39 

1.39 

16 

1.55 

1.42 

1.55 

1.54 

1.51 

1.51 

18 

1.43 

1.52 

1.57 

1.57 

1.58 

1.59 

20 

1.44 

1.54 

1.30 

1.30 

1.36 

1.36 
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Table 3 (Continued) 


X 

T 4 (-0.15) 

SCATTERING EFFICIENCY Q sca 
T 4 (0.15) T 6 (-0.05) T 6 (0.05) 

T 6 (-0.10) 

T 6 (0.10) 

2 

1.62 

1.61 

1.65 

1.65 

1.63 

1.63 

3 

3.02 

3.02 

3.06 

3.06 

3.06 

3.06 

4 

3.59 

3.58 

3.55 

3.55 

3.57 

3.57 

5 

3.28 

3.22 

3.24 

3.24 

3.25 

3.23 

6 

2.40 

2.46 

2.35 

2.35 

2.37 

2.37 

7 

1.68 

1.76 

1.54 

1.54 

1.62 

1.61 

8 

1.60 

1.61 

1.38 

1.38 

1.46 

1.46 

9 



— 

— 

— 

1.71 

1.72 

10 

2.11 

2.14 

2.00 

1.99 

2.04 

2.00 

12 

1.93 


1.79 

1.80 

1.84 

1.86 

14 



1.37 

1.38 

1.50 

1.51 

15 



1.44 

1.44 

1.55 

1.55 

16 



1.59 

1.59 

1.65 

1.65 


X 

T 8 (-0.05) 

T 8 (0.05) 

T 8 (-0.10) 

T 8 (0.10) 

2 

1.65 

1.65 

1.62 

1.61 

3 

3.06 

3.05 

3.05 

3.05 

4 

3.56 

3.56 

3.59 

3.61 

5 

3.24 

3.25 

3.23 

3.26 

6 

2.35 

2.35 

2.37 

2.36 

7 

1.53 

1.53 

1.59 

1.58 

8 

1.36 

1.37 

1.42 

1.43 


2.00 1.98 


2.00 
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Table 3 (Continued) 


X 

S 

Q 

avg 

ABSORPTION EFFICIENCY Q abs 
T 3 (0.05) T 3 (0.10) T 3 (0.15) 

T 2 (-0.10) 

T 2 (0. 10) 

2 

.170 

.170 

0.170 

0.169 

0.168 

0.169 

0.169 

3 

.259 

.261 

0.259 

0.259 

0.260 

0.259 

0.260 

4 

.352 

.357 

0.352 

0.352 

0.352 

0.355 

0.358 

5 

.471 

.455 

0.467 

0.458 

0.448 

0.456 

0.450 

6 

.518 

.515 

0.519 

0.518 

0.513 

0.514 

0.506 

7 

.543 

.573 

0.552 

0.566 

0.576 

0.565 

0.579 

8 

.641 

.628 

0.641 

0.642 

0.640 

0.625 

0.632 

9 

.684 

.664 

— 

— 

— 

0.679 

0.661 

10 

.695 

.702 

0.703 

0.713 

0.720 

0.701 

0.705 

12 

.729 

.740 

0.740 

0.764 

0.779 

0.745 

0.747 

14 

.765 

.769 

0.788 

0.810 


0.783 

0.782 

15 

.768 

.784 

0.795 

0.828 


0.793 

0.793 

16 

.808 

.800 

0.824 

0.845 


0.808 

0.809 

18 

.823 

.830 

0.848 

0.874 


0.841 

0.839 

20 

.865 

.861 

0.875 

0.897 


0.866 

0.866 


X 

T 2 (-0.20) 

T 2 (0.20) 

T 4 (-0.05) 

T 4 (0.05) 

T 4 (-0.1Q) 

T 4 (0.10) 

2 

0.167 

0.168 

0.170 

0.170 

0.169 

0.169 

3 

0.258 

0.264 

0.259 

0.259 

0.260 

0.260 

4 

0.357 

0.362 

0.353 

0.353 

0.355 

0.355 

5 

0.444 

0.434 

0.466 

0.466 

0.456 

0.455 

6 

0.519 

0.516 

0.516 

0.518 

0.509 

0.521 

7 

0.570 

0.574 

0.549 

0.551 

0.563 

0.563 

8 

0.623 

0.628 

0.641 

0.628 

0.644 

0.620 

9 

— 

— 

— 

— 

0.686 

0.697 

10 

0.708 

0.707 

0.695 

0.704 

0.711 

0.722 

12 

0.757 

0.763 

0.765 

0.761 

0.777 

0.776 

14 

0.797 

0.799 

0.782 

0.785 

0.824 

0.826 

15 

0.812 

0.814 

0.814 

0.807 

0.840 

0.842 

16 

0.829 

0.827 

0.831 

0.835 

0.854 

0.858 

18 

0.857 

0.853 

0.859 

0.856 

0.886 

0.892 

20 

0.881 

0.878 

0.887 

0.887 

0.913 

0.912 
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Table 3 (Continued) 


X 

T 4 (-0.15) 

ABSORPTION EFFICIENCY Q abs 
T 4 (0.15) T 6 (-0.05) T 6 (0.05) 

T 6 (-0.10) 

T 6 (0.10) 

2 

0.169 

mm 

0.170 

0.170 

0.169 

0.169 

3 

0.261 

WBm 

0.259 

0.259 

0.260 

0.260 

4 

0.356 


0.353 

0.353 

0.354 

0.354 

5 

0.449 

0.445 

0.468 

0.468 

0.462 

0.459 

6 

0.504 

0.526 

0.520 

0.521 

0.520 

0.526 

7 

0.572 

0.567 

0.549 

0.550 

0.564 

0.571 

8 

0.647 

0.636 

0.635 

0.632 

0.632 

0.629 

9 

— 

— 

— 

— 

0.697 

0.681 

10 

0.732 

0.738 

0.705 

0.704 

0.737 

0.735 

12 

0.789 


0.761 

0.761 

0.791 

0.777 

14 



0.795 

0.795 

0.835 

0.840 

15 



0.814 

0.816 

0.860 

0.857 

16 



0.835 

0.833 

0.879 

0.874 


X 

T g (-0.05) 

T 8 (0.05) 

Tg(-O.lO) 

T 8 (0.10) 

2 

0.169 

0.170 

0.168 

0.169 

3 

0.259 

0.259 

0.259 

0.259 

4 

0.353 

0.351 

0.350 

0.351 

5 

0.468 

0.467 

0.458 

0.459 

6 

0.522 

0.522 

0.526 

0.527 

7 

0.550 

0.550 

0.563 

0.566 

8 

0.631 

0.634 

0.624 

0.636 


9 — — 

0.702 


0.704 


0.711 
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Table 3 (Continued) 


SINGLE-SCATTERING ALBEDO 


S S avg T 3 (0.05) T 3 (0.10) T 3 (0.15) T 2 (-0.10) T 2 (0.I0) 


2 

.907 

.907 

0.907 

0.907 

0.907 

0.907 

0.907 

3 

.922 

.921 

0.922 

0.922 

0.921 

0.921 

0.921 

4 

.909 

.909 

0.910 

0.909 

0.909 

0.910 

0.909 

5 

.873 

.876 

0.874 

0.876 

0.878 

0.877 

0.879 

6 

.819 

.816 

0.819 

0.821 

0.825 

0.822 

0.823 

7 

.735 

.730 

0.734 

0.737 

0.750 

0.741 

0.737 

8 

.677 

.689 

0.683 

0.696 

0.716 

0.694 

0.695 

9 

.708 

.713 

— 

— 

— 

0.707 

0.713 

10 

.742 

.737 

0.738 

0.736 

0.738 

0.733 

0.729 

11 

.740 

.736 

— 

— 

— 

0.733 

0.731 

12 

.710 

.708 

0.706 

0.696 

0.687 

0.706 

0.709 

13 

.660 

.666 

— 

— ■ 


0.670 

0.673 

14 

.639 

.639 

0.631 

0.627 


0.645 

0.645 

15 

.647 

.645 

0.643 

0.638 


0.646 

0.643 

16 

.661 

.663 

0.654 

0.653 


0.656 

0.651 

17 

.670 

.667 

— 

— 


0.657 

0.657 

18 

.657 

.652 

0.648 

0.636 


0.645 

0.650 

19 

.625 

.625 

— 

— 


0.627 

0.633 

20 

.600 

.603 

0.596 

0.592 


0.613 

0.614 


X 

T 2 (-0.20) 

T 2 (0.20) 

T 4 (-0.Q5) 

T 4 (0.05) 

T 4 (-0.I0) 

T 4 (O.10) 

2 

0.906 

0.904 

0.907 

0.907 

0.907 

0.907 

3 

0.921 

0.919 

0.922 

0.922 

0.921 

0.921 

4 

0.910 

0.910 

0.910 

0.910 

0.910 

0.910 

5 

0.880 

0.884 

0.875 

0.874 

0.877 

0.877 

6 

0.827 

0.831 

0.820 

0.819 

0.823 

0.820 

7 

0.762 

0.765 

0.737 

0.738 

0.740 

0.746 

8 

0.716 

0.712 

0.683 

0.686 

0.695 

0.700 

9 

— 

— 

— 

— 

0.720 

0.715 

10 

0.721 

0.712 

0.742 

0.741 

0.742 

0.740 

11 

— 

— 

— 

— 

0.735 

0.738 

12 

0.699 

0.712 

0.702 

0.702 

0.705 

0.704 

13 

— 

— 

— 

— 

0.659 

0.655 

14 

0.668 

0.670 

0.631 

0.632 

0.623 

0.622 

15 

0.661 

0.645 

0.631 

0.633 

0.623 

0.623 

16 

0.652 

0.632 

0.651 

0.648 

0.639 

0.638 

17 

— 

— 

— 

— 

0.647 

0.646 

18 

0.625 

0.641 

0.646 

0.647 

0.641 

0.641 

19 

— 

— 

— 

— 

0.621 

0.621 

20 

0.620 

0.637 

0.594 

0.594 

0.598 

0 ; .599 
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Table 3 (Continued) 


X 

T 4 (— 0. 15) 

SINGLE-SCATTERING ALBEDO 
T 4 (0.15) T 6 (-0.05) T 6 (0.05) 

T 6 (-0.10) 

T 6 (0.10) 

2 

0.906 

0.906 

0.907 

0.907 

0.906 

0.906 

3 

0.920 

0.921 

0.922 

0.922 

0.922 

0.922 

4 

0.910 

0.910 

0.910 

0.910 

0.910 

0.910 

5 

0.880 

0.879 

0.874 

0.874 

0.876 

0.876 

6 

0.826 

0.824 

0.819 

0.819 

0.820 

0.818 

7 

0.746 

0.756 

0.737 

0.737 

0.742 

0.738 

8 

0.712 

0.717 

0.685 

0.686 

0.698 

0.699 

9 

0.729 

0.731 

— 

— 

0.710 

0.716 

10 

0.742 

0.744 

0.739 

0.739 

0.735 

0.731 

11 

0.737 

0.740 

— 

— 

0.733 

0.731 

12 

0.710 


0.702 

0.703 

0.699 

0.705 

13 

0.666 


— 

— 

0.664 

0.667 

14 



0.633 

0.634 

0.642 

0.643 

15 



0.639 

0.638 

0.643 

0.644 

16 



0.656 

0.656 

0.652 

0.654 

17 





0.650 

0.651 


X 

T g (-0.05) 

T 8 (0.05) 

Tg(-O.lO) 

T 8 (0.10) 

2 

0.907 

0.907 

0.906 

0.905 

3 

0.922 

0.922 

0.922 

0.922 

4 

0.910 

0.9.10 

0.911 

0.911 

5 

0.874 

0.874 

0.876 

0.877 

6 

0.818 

0.818 

0.818 

0.817 

7 

0.736 

0.736 

0.739 

0.736 

8 

n 

0.683 

0.684 

0.695 

0.692 

y 

10 

0.740 

0.740 


0.736 
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Table 3 (Continued) 


X 

S 

c 

^avg 

ASYMMETRY FACTOR g 
T 3 (0.05) T 3 (0.10) T 3 (0.15) 

T 2 (-0.10) 

T 2 (0.10) 

2 

.634 

.634 

0.634 

0.635 

0.636 

0.636 


3 

.747 

.747 

0.747 

0.746 

0.746 

0.746 

mm 

4 

.772 

.772 

0.771 

0.768 

0.766 

0.770 

1EG9 

5 

.751 

.752 

0.749 

0.745 

0.743 

0.751 


6 

.701 

.690 

0.700 

0.696 

0.691 

0.696 


7 

.611 

.604 

0.611 

0.613 

0.621 

0.618 

Wfim 

8 

.612 

.633 

0.614 

0.623 

0.640 

0.633 

0.634 

9 

.748 

.749 

— 

— 

— 

0.738 

0.739 

10 

.829 

.826 

0.822 

0.801 

0.786 

0.816 

0.815 

11 

.860 

.857 

— 

— 

— 

0.848 

0.849 

12 

.860 

.859 

0.854 

0.830 

0.801 

0.853 

0.857 

13 

.843 

.845 

— 

— 


0.845 

0.849 

14 

.838 

.839 

0.839 

0.820 


0.842 

0.844 

15 

.847 

.852 

0.855 

0.845 


0.854 

0.852 

16 

.868 

.871 

0.874 

0.867 


0.868 

0.866 

17 

.884 

.882 

— 

— 


0.878 

0.877 

18 

.884 

.883 

0.886 

0.878 


0.880 

0.882 

19 

.876 

.879 

— 

— 


0.880 

0.882 

20 

.875 

.878 

0.877 

0.874 


0.881 

0.882 


X 

T 2 (-0.20) 

T 2 (0.20) 

T 4 (-0.05) 

T 4 (0.05) 

T 4 (-0.10) 

T 4 (0.10) 

2 

0.638 

0.645 

0.636 

0.636 

0.639 

0.639 

3 

0.741 

0.752 

0.749 

0.749 

0.751 

0.752 

4 

0.771 

0.775 

0.772 

0.772 

0.774 

0.775 

5 

0.757 

0.762 

0.751 

0.750 

0.750 

0.750 

6 

0.707 

0.716 

0.696 

0.696 

0.690 

0.688 

7 

0.653 

0.666 

0.606 

0.607 

0.600 

0.611 

8 

0.659 

0.657 

0.613 

0.613 

0.617 

0.618 

9 

— 

— 

— 

— 

0.727 

0.719 

10 

0.792 

0.780 

0.820 

0.817 

0.804 

0.798 

11 

— 

— 

— 

— 

0.835 

0.832 

12 

0.834 

0.851 

0.852 

0.852 

0.834 

0.833 

13 

— 

— 

— 

— 

0.818 

0.813 

14 

0.851 

0.859 

0.834 

0.834 

0.811 

0.805 

15 

0.861 

0.857 

0.850 

0.850 

0.829 

0.826 

16 

0.868 

0.860 

0.875 

0.874 

0.856 

0.852 

17 

— 

— 

— 

— 

0.874 

0.870 

18 

0.874 

0.882 

0.890 

0.890 

0.881 

0.880 

19 

— 

— 

— 

— 

0.882 

0.881 

20 

0.889 

0.893 

0.886 

0.885 

0.881 

0.880 
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Table 3 (Continued) 


X 

T 4 (-0.15) 

ASYMMETRY FACTOR g 
T 4 (0.15) T 6 (-0.05) T 6 (0.05) 

T 6 (-0.10) 

T 6 (0.10) 

2 

0.644 

0.645 

0.637 

0.637 

0.644 

0.644 

3 

0.754 

0.755 

0.749 

0.749 

0.753 

0.753 

4 

0.775 

0.778 

0.773 

0.773 

0.777 

0.777 

5 

0.751 

0.751 

0.752 

0.753 

0.754 

0.753 

6 

0.688 

0.687 

0.698 

0.698 

0.691 

0.693 

7 

0.603 

0.622 

0.607 

0.607 

0.602 

0.607 

8 

0.626 

0.630 

0.605 

0.607 

0.605 

0.614 

9 

— 

— 

— 

— 

0.705 

0.712 

10 

0.789 

0.781 

0.817 

0.817 

0.788 

0.787 

1! 

0.818 

0.813 

— 

— 

0.819 

0.821 

12 

0.820 


0.845 

0.845 

0.809 

0.826 

13 

0.808 


— 

— 

0.816 

0.816 

14 



0.838 

0.836 

0.822 

0.820 

15 



0.853 

0.852 

0.841 

0.840 

16 



0.875 

0.873 

0.862 

0.860 


X 

T g (-0. Q 5) 

T 8 (0.05) 

T 8 (-0.10) 

T 8 (0.10) 

2 

0.637 

0.637 

0.644 

0.647 

3 

0.750 

0.750 

0.756 

0.756 

4 

0.774 

0.774 

0.776 

0.778 

5 

0.751 

0.751 

0.750 

0.753 

6 

0.700 

0.701 

0.694 

0.694 

7 

0.613 

0.614 

0.607 

0.608 

8 

0.610 

0.610 

0.602 

0.603 


0.818 0.787 


10 


0.820 



X. SCATTERED INTENSITIES IN FIXED ORIENTATION 


One goal of this research was to provide a set of results against which the validity of the EBCM could 
be tested, either by other exact theoretical techniques or by laboratory experiments. All our random- 
orientation results are available for this purpose. However, they may not be the best choice, for 
several reasons: 

® random-orientation intensities are rather smooth functions of angle; correctness of the 
EBCM could be more incisively tested against more highly structured angular patterns, dif- 
fering more from a sphere; 

• random-orientation results depend to some extent on the quadrature over orientation, which 
introduces an unknown error; 

• comparison with microwave-analogue measurements (see Appendix A) would be easier in 
fixed orientation; 

• it is intensities which are measured in the laboratory, while all our random-orientation plots 
are for phase function and degree of polarization; it requires some effort to deconvolve in- 
tensities from these plots 

® our random-orientation results are for unpolarized incident radiation, whereas laboratory- 
sources (esp. lasers) are frequently polarized. 

Hence, we have decided to include a small selection of fixed-orientation intensities, using both parallel 
and perpendicular polarized incident light, in Appendix B. These are for Chebyshev particles with e = 
.0.10 and x = 5,10. Shapes were chosen so that the EBCM could be tested for convex |[T 2 ], concave 
[T 3 , T 4 , T 6 ], and asymmetric [T 3 ] particles. Two fixed orientations for which no cross-polarization oc- 
curs were chosen: nose-on and ‘perpendicular-to-nose-on’ (9 p = 90°, 4> p = 0). 

Solid lines in App. B refer to the equal-volume sphere; dotted lines to the nose-on orientation; and 
dashed lines to the ‘perpendicular-to-nose-on’ orientation. The ‘fair’ convergence criterion was not 
used to produce any of these results. 
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XI. PHASE FUNCTION AND DEGREE OF POLARIZATION IN RANDOM ORIENTATION 


MIXTURE 

Figures 9 and 10 show the phase function (55) and degree of polarization (58), respectively, for the 
MIXTURE of 23 particles. These are the most general results of our study, vis a vis the ‘net’ effect of 
nonsphericity plus random orientation on the angular variation of scattering. 

The arrangement of these plots is similar to that in Appendices C and D (described below). Each phase 
function plot (for 60-180°) is on the same page with its corresponding percentage-difference-from-a- 
sphere plot (for 0-180°). The MIXTURE phase function for x = 2,3,4,5,6,8,10,15, and 20 (solid line) 
can be compared with the exact (dotted line) and Ax = 0.1x size-averaged (dashed line) spherical phase 
functions. The solid, dashed, and dotted lines have the same meaning on the degree of polarization 
plots (Fig. 10). The percentage differences in Fig. 9 are with respect to the exact (dotted) and size- 
averaged (dashed) Mie phase functions. 

Although we calculated results for many more integer values of ‘x’ in the range 10-25 than we show 
here, the decision was made to omit these results because, too often, the changes from one value of ‘x’ 
to the next were not noteworthy. 

The two spherical curves in Figure 9 bear out our remarks in Sec. VIII, to the effect that down-spikes 
in the angular range 80-150° tend to be washed out by size-averaging. One should be cautious of the 
visual impression that size-averaging biases the phase function upward, however; often, this is due on- 
ly to plotting on a logarithmic vertical scale. 

T n particles 

The bulk of our results for phase function (60-180° only) and degree of polarization are contained in 
Appendices C and D, respectively. On each such plot, the solid line shows the Ax = 0.1x size-averaged 
spherical result, while the dotted, dashed, and dot-dash lines refer to T n particles with successively in- 
creasing values of |e|. Note the change of convention relative to the MIXTURE; there, we used a 
solid line for the nonspherical result because we wanted to show two spherical results (exact and size- 
averaged). 

These Appendices were designed to allow the reader to flip through them rapidly, generating almost a 
moving picture of spherical-nonspherical differences as a function of particle size. Thus, a rapid feel- 
ing for the results can be obtained. In Appendix C, the phase functions are plotted at the top of each 
page, and their percent differences from a sphere at the bottom; thus, the ‘action’ can be stopped at 
any point to compare the two. 
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MIXTURE x=2. 



MIXTURE x = 2. 
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Figure 9. 60-180 degree phase function for nonspherical MIXTURE (solid line), exact sphere 
(dotted line), and Ax = O.lx size distribution of spheres (dashed line) for x = 
2,3,4,5,6,8,10,15, and 20. Companion plots show percent differences of the' 
nonspherical phase functions from the spherical ones for the full range of angles 
(0-180 degrees); interrupted lines refer to the same spherical cases as in the 
corresponding phase function plot. Note how size-averaging the spherical results 
frequently mimics the effect of nonsphericity. 
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MIXTURE x=3. 
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MIXTURE x=3. 



Figure 9 (Continued) 
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MIXTURE x=5. 



Figure 9 (Continued) 
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MIXTURE x=6. 
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MIXTURE x=8. 
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Figure 9 (Continued) 
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MIXTURE x= 15. 



MIXTURE x= 15. 



Figure 9 (Continued) 
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MIXTURE x = 20. 
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MIXTURE x = 20. 



Figure 9 (Continued) 
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SCATTERED INTENSITIES IN FIXED ORIENTATION 
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MIXTURE x = 2. 



MIXTURE x = 3. 



SCATTERING ANGLE 

Figure 10. Degree of polarization for nonspherical mixture (solid line), exact sphere (dotted line), 
and Ax = O.lx size distribution of spheres (dashed line) for x= 2,3,4,5,6,8,10,15 
and 20. Note how, as in Fig. 9, size-averaging the spherical results frequently mimics 
the effect of nonsphericity. 
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MIXTURE x=5. 



Figure 10 (Continued) 





DEGREE OF POLARIZATION OEGREE OF POLARIZATION 


SCATTERED INTENSITIES IN FIXED ORIENTATION 


87 



0 20 40 60 80 100 120 140 160 190 

SCATTERING ANGLE 


1 .0 
.8 
.6 
.4 
.2 
0 

-.2 
-.4 
- .6 
- .0 
1 .0 

o 20 40 60 00 100 120 140 160 180 

SCATTERING ANGLE 



Figure 10 (Continued) 
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MIXTURE x=10. 



MIXTURE x=15. 



Figure 10 (Continued) 
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MIXTURE x = 20. 
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Figure 10 (Continued) 
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The organization of the Appendices can be likened to a sequence of 3 nested loops. The idea was, to 
group the most similar results closest together, then the next most similar, and so forth. Thus, the in- 
nermost ‘loop’ is over the magnitude of the deformation parameter e of a Chebyshev particle; results 
from those values of e for which the EBCM converged are plotted on the same frame. 

The reason for grouping different values of e closest together, is to test our intuitive notions of pertur- 
bation theory. When e doubles, so these notions go, the differences from the spherical phase function 
and degree of polarization should double also. Thus, dividing the percent differences by the cor- 
responding | e | should cause them to fall almost on top of one another. This does not in fact happen 
— far from it — else we should have presented the percentages in this way. 

The middle ‘loop’ is over particle type. Since T n + and T n _ particles bear a subtle similarity, as describ- 
ed in Sec. II, they are placed on facing pages, for easy comparison. The first particle in each group is 
T 3 , since it is the only unsymmetric one; then the other T n ’s follow in order of n 2, 4, 6, 8. T 20 
results are omitted because only the | e | = 0.05, x < 8 cases could be converged, and they did not dif- 
fer much from a sphere. 

The outermost loop is over size parameter x (from 2 to 20), because the biggest changes occur as size 
changes. As noted above, we omitted all of our results from x = 11 to x = 25 except x = 15 and 20, 
because we were unable to see big enough differences among the x > 10 plots — just an increasingly 
complex jumble of oscillations. We feel that further filtering — perhaps the subtraction of a common 
oscillatory term from all curves — is necessary to make any sense out of this data. Or perhaps phase 
function and degree of polarization are the wrong quantities to look at; they just do not make 
spherical-nonspherical differences stand out in stark relief. 

Until a better way is found to look at spherical-nonspherical differences for x > 10, our use of the 
EBCM all the way up to x = 25 (no one had gone beyond x= 10 before) remains more of a tours de 
force than a way to gain fundamentally new insights. 

We tried to keep the vertical-axis range of our plots as small as possible, in order to magnify spherical- 
nonspherical differences. At the same time, we wanted to avoid having a different scale for every plot. 
The compromise we struck was to fix the range for large groups of plots (for example, at one full 
decade), but to slide this range up or down to match the quantities being plotted, as necessary. This 
way, the spread between curves on neighboring plots is comparable, even when the vertical scales are 
not identical. 

It was also the question of vertical resolution which drove us to plot phase function separately for 
0-60° and 60-180°. The forward peak of the phase function takes several decades to contain it, yet 
spherical-nonspherical differences are all but negligible there. Meanwhile, the important differences 
in side-scattering were shrunk to Lilliputian proportions. 

Note that the percentage difference plots in App. C change from a linear to a logarithmic vertical scale 
past x = 6. This is because the percentages become very large — over 100% in many cases — in spite 
of the relatively small deviations from sphericity. Percentages below about 1%, we regard as in the 
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noise level, and have explicitly accounted for this in the logarithmic plots by starting the scale at 1% 
(which is labelled ‘0’). 

0-00° phase f unction 

We present in Fig. 1 1 only a sparse sampling of our phase function results in the angular range 0-60°. 
(The corresponding percentage differences from a sphere for 0-60° can be found in Appendix C.) The 
meanings of the various lines are the same as in Fig. 9 (for the MIXTURE cases) and as in Appendix C 
(for the T n cases). 

The reason for omitting much of this data, is simply that it proved uninstructive. After all, 0-60° is 
primarily a diffraction region, dominated by ‘rays’ that do not actually ‘hit’ a particle but merely bend 
around it. (One cannot really talk of ‘rays’ at x = 1 , but already by x = 10 the concept has some utility.) 
These ‘rays’ are insensitive to particle shape and refractive index. They are sensitive primarily to pro- 
jected area (van de Hulst, 1981); and the differences in projected area between randomly-oriented 
Chebyshev particles and equal-volume spheres are minute. 

Hence, as might be expected, the spherical and nonspherical phase functions tracked each other very 
closely in this region — so closely, often, that they were indistinguishable. Some of the plots in Fig. 1 1 
were selected to illustrate this point. The remainder show cases where the differences are more 
substantial. 

Acknowledgments: We are indebted to the National Center for Atmospheric Research for providing 
the computer time for these studies; to David Kennison of NCAR for assistance with some of the 
graphical presentations; and to NASA for publishing this compendium, whose completion spanned 
over 5 years and hence came to be affectionately known as ‘The Eternal Project’. 
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Figure 11. 0-60 degree phase function for nonspheres (solid Sine for MIXTURE, interrupted lines 
in other cases — dotted, dashed, and dot-dash lines correspond to the various 
e values shown, in increasing order) and for Ax = 0.1 x size distribution of spheres 
(dashed line for MIXTURE, solid line in other cases). Exact sphere result (dotted line) is 
shown only for MIXTURE. Size parameters x = 5, 10, 15, and 20 o 
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(V) le!=0. 10,0.20 x=5. 



(V) Id — 0 . 1 0,0.20 x = 5. 



Figure 1 1 (Continued) 
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Figure 1 1 (Continued) 
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(v) l£l=0.05,0. 10 x=5. 



(V) lcl=0.05,0. 10 x-5. 



Figure 1 1 (Continued) 
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Figure 1 1 (Continued) 
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Figure 11 (Continued) 
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MIXTURE x= 15. 
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Figure 11 (Continued) 
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(V) lel=0. 10,0.20 x= 15. 



(V) lel=0. 10,0.20 x- 15. 



Figure 11 (Continued) 
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(V) lel=0.05,0.10 x= 15. 



(T e + > lcl=0.05,0. 10 x= 15. 



Figure 11 (Continued) 
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(T a -> |£|=0. 10,0.20 x=2Q. 



<T 2 + ) lel=0. 10,0.20 x=20. 



Figure 11 (Continued) 
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Figure 11 (Continued) 





REFERENCES 


Asano, S. and G. Yamamoto, 1975: Light scattering by a spheroidal particle, Appl. Opt, 14 , 29 - 49 , 

Asano, S., 1980: Light scattering by randomly oriented spheroidal particles, Appl. Opt. 19 , 962 - 974 . 

Barber, P., 1973: “Differential scattering of electromagnetic waves by homogeneous isotropic dielec- 
tric bodies”, Ph.D. Thesis, UCLA, Los Angeles, Calif, (available from University Microfilms, Ann 
Arbor, Mich.) 

Barber, P., 1977: Resonance electromagnetic absorption by nonspherical dielectric objects, IEEE 
Transac. Microwave Theory & Tech. MTT-25, 373-381. 

Barber, P. and C. Yeh, 1975: Scattering of electromagnetic waves by arbitrarily shaped dielectric 
bodies, Appl. Opt. 14, 2864-2872. 

Bowman, J., Senior, T. and P. Uslenghi, eds., 1969: Electromagnetic and Acoustic Scattering By Sim- 
ple Shapes, American Elsevier, New York. 

Chylek, P., V. Ramaswamy and W. Wiscombe, 1982: Note on “The scattering of radiation by 
moderately nonspherical particles”, I. Atmos. Sci. 39, 1886-1889. 

Collis, R. and P. Russell, 1976: Lidar measurement of particles and gases by elastic backscattering and 
differential absorption, in Laser Monitoring of the Atmosphere , E. Hinkley, ed., Springer-Verlag, 
New York, pp. 71-151. 

Dave, I., 1969: Effect of coarseness of the integration increment on the calculation of the radiation 
scattered by polydispersed aerosols, Appl. Opt. 8, 1161-1167. 

Davis, P. and P. Rabinowitz, 1975: Methods of Numerical Integration, Academic Press, New York, 
258 pp. 

Gerber, H., 1979: Absorption of 632.8 nm radiation by maritime aerosols near Europe, J. Atmos. Sci. 
36, 2502-2506. 

Joseph, I., W. Wiscombe and J. Weinman, 1976: The delta-Eddington approximation for radiative 
flux transfer, J. Atmos. Sci. 33, 2452-2459. 

Lentz, W., 1975: Generating Bessel functions in Mie scattering calculations using continued fractions, 
Appl. Opt. 15, 668-671. 


109 


PREQEDO^Q PAGE BLANK NOT FILMED 



110 


SINGLE SCATTERING FROM NONSPHERICAL CHEBYSHEV PARTICLES 


Millar, R., 1969: Rayleigh hypothesis in scattering problems, Elec. Lett. 5, 416-418. 

Millar, K., 1973: The Rayleigh hypothesis and a related least-squares solution to scattering problems 
for periodic surfaces and other scatterers, Radio Sci. 8, 785-796. 

Morse, P. and H. Feshbach, 1953: Methods of Theoretical Physics, McGraw-Hill, New York. 

Mugnai, A. and W. Wiscombe, 1980: Scattering of radiation by moderately non-spherical particles, 
J. Atmos. Sci. 37, 1291-1307. 

Nussenzveig, H. and W. Wiscombe, 1980: Forward optical glory, Opt. Lett. 5, 455-457. 

Patterson, T., 1973: Algorithm for automatic numerical integration over a finite interval, Commun. 

ACM 16, 694-699. 

Pollack, I. and J. Cuzzi, 1980: Scattering by nonspherical particles of size comparable to a 
wavelength: a new semi-empirical theory and its application to tropospheric aerosols, J. Atmos. Sci. 

37, 868-881. 

Pruppacher, H. and R. Fitter, 1971: Calculations and measurements of the shape of deformed rain- 
drops, 3. Atmos. Sci. 28, 86-94. 

Ross, W., 1972: Computation of Bessel functions in light scattering studies, Appl. Opt. 11, 

1919—1923. 


Schelkunoff, S., 1943: Electromagnetic Waves , van Nostrand, New York. 

Schuerman, D., ed., 1980: Light Scattering by Irregularly Shaped Particles, Plenum Press, New York. 

Schuerman, D., R. Wang, B. Gustafson and R. Schaefer, 1981: Systematic studies of light scattering. 

1: Particle shape, Appl. Opt. 20, 4039-4050. 

Shipley, S. and J. Weinman, 1979: A numerical study of scattering by large dielectric spheres, J. Opt. 

Soc. Amer. 68, 130-134. 

Van de Hulst, H., 1981: Light Scattering by Small Particles, Dover, New York. 

Waterman, P., 1965: Matrix formulation of electromagnetic scattering, Proc. IEEE 53, 805-812. 

Waterman, P., 1971: Symmetry, unitarity, and geometry in electromagnetic scattering, Phys. Rev. 

D3, 825-839. 

Wiscombe, W., 1979: Mie Scattering Calculations: Improvements in Technique and Fast, Vec- 
tor-Speed Computer Codes, NCAR Technical Note TN-140 + STR, National Center for At- 
mospheric Research, Boulder, Colorado. 



REFERENCES 


Wiscombe, W., 1980: Improved Mie scattering algorithms, Appl. Opt. 19, 1505-1509. 

Wiscombe, W. and A. Mugnai, 1980: Exact calculations of scattering from moderately-nonspherical 
Chebyshev particles: comparisons with equivalent spheres, in Light Scattering by Irregularly Shaped 
Particles , D. Schuerman, ed., Plenum Press, New York, pp. 141-152. 

Zerull, R., 1976: Scattering measurements of dielectric and absorbing nonspherical particles, Beitr. 
Phys. Atmos. 49, 166-188. 

Zerull, R., Giese, R. and K. Weiss, 1977: Scattering functions of nonspherical dielectric and absorbing 
particles vs. Mie theory, Appl. Opt. 16, 777-778. 




APPENDIX A: REVIEW OF NONSPHERICAL SCATTERING 


There is barely a field of science or engineering that does not have some interest in scattering of radia- 
tion. In many of these fields, the problem can be separated into single scattering by individual objects, 
followed by multiple scattering among them. But rarely do scattering objects have those ideal shapes 
— sphere, spheroid, circular cylinder, and so forth — that allow an ‘exact’ solution of Maxwell’s 
Equations in infinite series of eigenfunctions (Morse and Feshbach, 1953, pp. 492-523). 

Instead, the problem of scattering by rather nasty-looking shapes constantly presents itself. Mostly, 
the effect of shape is ignored; as van de Hulst (1981) says, “the formulae for spherical particles are us- 
ed in 95 Vo of all applications.” Assuming sphericity is certainly the path of least resistance; 
well-documented computer algorithms for spherical (‘Mie’) scattering are available (Wiscombe, 1979, 
1980) and have been available since Dave produced the first ones in 1968. Circular cylinder and 
spheroidal algorithms are also available, although less widely used. 

But the assumption of sphericity is rarely made after first studying the effect of nonsphericity, then 
discounting it. Instead, it is almost invariably made a priori. Hence, there is an increasing desire to go- 
back and check this assumption, and find out where it breaks down. 

The questions one asks about nonspherical scattering depend, of course, on one’s research interests. 
Those who use lidar and radar to sound the atmosphere, for example, want to know only about the in- 
tensity and polarization in direct backscatter . Military planners want to block enemy radar and com- 
munications with airborne ‘chaff’, whose shape (among other factors) is adjusted to obtain the max- 
imum extinction cross-section with the minimum amount of material. 

Medical researchers want to side-scatter laser light off blood samples and determine the state of the 
suspended cells. And they want to know the absorption cross-section of human organs in the 
microwave region. Aerosol specialists want to infer aerosol size distributions from the angular 
distribution of near-forward-scattered light. Climate theorists want to know how ice crystals in cirrus 
clouds control the Earth’s radiative fluxes . 

Each group has a different set of questions it wants answered. We shall begin by reviewing the ques- 
tions that are being asked about nonspherical scattering, in order to set the stage for the following 
discussions of experimental and theoretical results. 

A.1 CURRENT QUESTIONS 

Neglecting quantum and non-linear effects, the questions revolve around one or more of the follow- 
ing list of classical scattering quantities, arranged in order of increasing complexity: 

extinction cross-section 

scattering cross-section 
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absorption cross-section 
single-scattering albedo 
radiation pressure cross-section 
asymmetry factor 

unpolarized scattered intensity at any angle from the direction of the incident radiation 

unpolarized phase function 

Legendre moments of unpolarized phase function 

polarized scattered intensities, parallel and perpendicular to the plane of scattering (the plane of 
the incident and scattered beams) 

Stokes parameters 

Mueller matrix 

Legendre moments of Mueller matrix elements 

To add to this list, there are different possible states of the incident radiation: plane wave, spherical 
wave, finite beam, and so forth. And while the above quantities all traditionally refer to the ‘far-field’ 
(far enough from the particle that radial components of the scattered field can be neglected), there are 
also important problems of near-field scattering as well — for example in paint or snow. 

The orientation of the particle is another complexifying factor. Some particles are simply in fixed 
orientations, as in microwave-analogue experiments. Others are partially randomly oriented, as ice 
needles in a cirrus cloud, or falling raindrops. Most are completely randomly oriented, as aerosol 

particles. 

Generally, as one moves down the above list, the effect of nonsphericity may be expected to become 
more and more important. Those who are interested only in cross-sections may be justified in ignoring 
nonsphericity in 95% of the cases. On the other hand, those who need the lower right-hand 2x2 sub- 
matrix of the Mueller matrix may find it almost impossible to ignore nonsphericity. Those whose 
problem areas lie in between, are in a grey area where they may or may not ignore it, depending on the 
size, deviation from sphericity, and refractive index of their scattering objects. 

Refractive index plays a particularly important role. When the real index differs significantly from 
unity, and the imaginary index is nearly zero, nonsphericity exerts its maximum effect. On the other 
hand, when real index is near unity (as for a bacterium suspended in water), or the imaginary index is 
large (as for a carbon particle), nonsphericity can often be safely neglected. 
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Questions of nonsphericity may also be embedded in a larger context: for example, multiple scatter- 
ing, or size-averaging. Size-averaging and shape-averaging wash out the fine details of nonspherical 
scattering, making the ensemble (we assume) look more like a collection of spheres. And the greater 
the extent of multiple scattering, the less the fine details of the phase function matter; only the cross- 
sections and the first few Legendre moments play a role. In such settings, nonsphericity may be more 
or less important, depending on the width of the size distribution and the total optical thickness of the 
scattering medium. 

Looking at the list of scattering quantities above, and the variety of complexifying factors, it is clear 
that a program of research in nomspherical scattering must be carefully focussed to address very 
specific questions. Otherwise, it risks being swamped by the sheer number of possible variables, and 
the variety of parameters which can influence those variables. Up to the present, the foci have been: 

• Can the nonspherical particle be replaced by an ‘equivalent sphere’? If so, which equivalent 
sphere (equal-volume, equal-area, . . .) is best? 

• Can we make simple, empirical adjustments to Mie theory to mimic nonspherical scattering? 

• When Mie theory just won’t work, can we use either an exact solution for cylinders or 
spheroids, or a small- or large-particle approximation? 

• Can we discern general features of nonspherical scattering just by looking at a few canonical! 
shapes (in the lab or theoretically)? 

® What impact does nonsphericity have on our remote sensing capabilities, ranging from in- 
terstellar dust effects on astronomical observations to ice particle effects on radar? 

The general tenor of the questions in nonspherical scattering is therefore somewhat as follows: What is 
the least amount of shape information we need to know? Can we just make simple modifications to 
existing theories to get what we want? For the fact of the matter is, few people outside of elecljical 
engineering (where unusual antenna shapes are routinely studied) want to know the exact scattering 
from a doughnut or a dodecahedron in fixed orientation. They just want to know the genera! effect of 
nonsphericity, shorn of its details. In what follows, this is the tack we shall take in describing ex- 
periments and theoretical results. 

A.2 MEASUREMENTS 

We will survey only a selection of measurements specifically aimed at studying nonsphericity. There 
are a large number of instruments which routinely measure the scattering from nonspherical particles. 
Yet these are of no use for our purposes, because they all make simple empirical adjustments to ac- 
count for nonsphericity (if they account for it at all). Thus, their goal is not to study nonsphericity, but 
to get rid of it, like an unwelcome dinner guest. 
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Traditionally, the measurements have tended to fall into two categories: 

• visible light scattering from micron-sized particles 

• microwave scattering by centimeter-sized objects 

In the future, much more of the electromagnetic spectrum will be available, using tuneable lasers and 
cryogenically cooled detectors; but until recently, experimenters were strongly constrained by the 
relatively primitive state of source and detector technology, and by the lack of windows in the spec- 
trum of the Earth’s atmosphere. 

Visible light scattering experiments involve relatively simple and cheap instrumentation, and have just 
as frequently been done in the field as in the laboratory. Microwave scattering experiments, by con- 
trast, require expensive instrumentation and large laboratories, and have as a consequence been 

undertaken by only a few groups. 

Microwave 

The “microwave-analogue method” was pioneered by Greenberg (1960, 1961), with the application 
of interstellar dust in mind. The idea was, to manufacture centimeter-sized scattering objects with 
desired shapes and refractive indices. Then, in a large, anechoic chamber, perhaps 100 feet across, one 
directed a microwave beam at the object and measured the scattering at all angles. By special tech- 
niques, even the extinction cross-section could be measured; and for independent confirmation, one 
could place temperature sensors in the object and measure the Joule heating directly (to get absorption 
cross-section). 

The ‘analogue’ part of the method meant simply that the results could be extrapolated to smaller or 
larger sized scatterers, keeping the ratio size/wavelength fixed, in much the same way as 
aerodynamicists extrapolate wind tunnel model results to full-sized aircraft. 

Very few published results came out of Greenberg’s effort. Greenberg, et. al. (1971) presented a few 
measurements of extinction cross-section for spheroids, while Greenberg (1972) made measurements 
for stacked cylinder configurations. Wang (1980) has given extinction cross-sections for aggregates of 
2, 4, and 8 spheres, as well as stacked cylinders. 

Waterman (1965, 1971) measured scattered microwave intensities from hemispherically capped 
cylinders and cones made of metal, primarily in order to verify his EBCM calculations. 

Zerull (1976) and Schuerman, et. al. (1981) have published by far the largest collection of microwave 
analogue measurements. They measured scattered intensities from roughened spheres, cubes, oc- 
tahedra, and irregular convex and concave particles, with various effective size parameters up to 20. 
Comparing to equal-volume spheres, their main conclusions were that: 

® for size parameters x < 6, agreement is within a factor of two; as x increases beyond 6, dif- 
ferences escalate dramatically 
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• nonspherical curves have greatly damped oscillations as a function of angle 

® agreement is very good in the forward peak (scattering angles 0-30 °), with nonspherical par- 
ticles probably scattering somewhat less than spheres there 

® nonspherical particles side-scatter (angles 30-140°) as much as an order of magnitude more 

• transparent nonspherical particles backscatter less , and their backscattering exhibits a much 
milder variation with angle 

• backscatter from opaque nonspherical particles shows an increase toward 180°, while the 
backscattering for spheres shows no such rise 

Zerull also found that concave particles tended to scatter more energy at most angles than 
equal-volume convex particles. 

Visible 

Visible light measurements have traditionally suffered from an inability to calibrate the instrumenta- 
tion using particles of known shape, size and refractive index. This problem has been somewhat 
alleviated of late, however, as monodisperse polystrene latex spheres have been more reliably 
manufactured (e.g.. Bottiger, et. ah, 1980). Also, the ‘vibrating-orifice aerosol generator 5 can 
generate quasi-monodispersions of various water-soluble aerosols (Pinnick, et. ah, 1976; Coletti, 
1984). 

There are many more visible than microwave measurements, reflecting the relative simplicity and com- 
pactness of sources and detectors at visible wavelengths; but these measurements are almost exclusive- 
ly of scattered intensity in the 10-170 degree angular range. The arrangement of source and detector 
precludes measurements near the forward or backward directions, although Ashkin and Dziedzic 
(1980, 1981), in their ‘optical levitation’ technique involving suspension of particles by laser light 
pressure alone, have obtained direct backscatter measurements. 

Lacking especially the important 0-10 degree scattering, it is difficult to obtain an accurate estimate of 
scattering cross-section (and hence of phase function) from integrating the scattered intensities over 
angle. As a result, there has been a tendency to normalize experimental results to the equal-volume 
sphere phase function at 10°, which may introduce biases into some of the conclusions drawn. Recent- 
ly, however, Coletti (1984) has used Fraunhofer diffraction theory to extrapolate his results into 
0-10°, which goes far toward correcting this problem. 

Extinction cross-section measurements have suffered from the classic problem that a detector with a 
finite aperture picks up some of the diffracted light. Depending on how much of this light is detected, 
the extinction can be mis-estimated by as much as a factor of two. Careful correction for the diffrac- 
tion can eliminate this effect if the particle projected area is known, but too often this is not the case. 
And with large errors in both extinction and scattering cross-section, almost nothing can be said about 
absorption cross-section, and hence the single-scattering albedo. 
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A number of early results on nonspherical scattering were reported at the First International Con- 
ference on Electromagnetic Scattering in 1963. Hodkinson (1963) measured extinction and for- 
ward-hemisphere scattering for quartz, diamond, flint, and coal dust particles of sizes 1 to 6 microns, 
in aqueous suspension. The curve of extinction as a function of size parameter x (= 10-50) did not 
display the large oscillations found in its Mie theory counterpart; instead, it just rose monotonically to 
an asymptote. Fraunhofer diffraction plus geometrical optics for a sphere gave a fair fit to the results 
for 0-60°, but not for 60-90°. 

Napper and Ottewill (1963) measured scattered intensities from transparent cubes and octahedra of 
size 0.5-0. 7 microns. The octahedral and equal-volume Mie results agreed well in 60-120 degrees, but 
the cubes side-scattered much more in both polarizations. 

Huffman and Thursby (1969) measured scattered intensities from irregular ice crystals, as well as hex- 
agonal plates and columns, in the angular range 10-150°. They showed that if they picked a size of 
sphere with equal scattering at 10 degrees, that sphere scattered considerably less to the side (around 

90°). 

In a well-done classic study, Holland and Gagne (1970) measured the Stokes parameters of light scat- 
tered from 0.1 to 1 micron plate-like, randomly oriented quartz crystals. They argue for comparing to 
equai-projected-area rather than equal-volume spheres, which gives good agreement out to 40° but 
underestimates the nonspherical side-scattering from there to 140 degrees or so. Past 140°, the situa- 
tion was reversed: Mie theory was more than a factor of 3 higher. Waggoner, et. al. (1972) also found 
that notisphericity depressed backscattering, although it must be emphasized that neither study could 
reach 180°. 

Proctor and Harris (1974) and Proctor and Barker (1974) measured extinction for numerous size 
distributions (spanning 0.1 to 50 microns) of quartz and diamond dust suspended in water. When they 
plotted extinction vs. the phase shift parameter 

p = 2 x (m — 1) 

they found no oscillations, unlike Mie theory; only the first main peak showed up. 

Chylek, et. al. (1976) measured scattered intensities for transparent salt particles of size 0.04 to 2 
microns. They showed a result which seemed contrary to everyone else’s, namely less side-scattering 
for non-spheres than for spheres. However, it seems that they used polarized incident light, whereas 
the previous measurements used unpolarized light. In fact, for one of the two incident polarizations in 
Holland and Gagne’s experiment, a similar result to Chylek, et. al. was obtained. Pinnick, et. al. 
(1976) continued this work, and found that either assuming voids in the particles, or averaging over 
sizes in Mie theory, improved agreement with spherical results. 

Perry, et. al. (1978) made one of the rare studies of the complete Mueller matrix. They used cubic salt 
and rounded ammonium sulfate particles of size 0.1 to 2 microns. They found that Mie theory could 
fit results for the rounded particles all the way up to x = 12, but for the cubes only up to x = 3. The 
nonspherical particles always had higher side-scatter, and lower backscatter, than equal-volume 

spheres. 
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Bottiger, et. al. (1980) also measured the Mueller matrix. Their scatterers were not ‘natural’ aerosol 
particles, but rather aggregates of 2, 3, or 4 polystyrene latex spheres, levitated electrostatically. They 
thus achieved a much greater degree of shape control than in previous experiments. They found that 
the strong angular oscillations in spherical Mueller matrix elements tend to be washed out — the more 
so, the more spheres were in the nonspherical aggregate. 

Saunders (1980) suspended 50 micron salt crystals on spider webs and measured scattered intensity at 
180° as they hydrated. He also subjected the hydrating particle to an electric field, causing distortions 
undetectable by microscope yet sufficient to cause large variations in backscattered intensity. This 
proved that direct backscatter is very sensitive to even slight nonsphericity. 

Coletti (1984) has scattered laser light from quasi-monodispersions of transparent salt particles and 
absorbing dye particles. He finds almost a complete washing-out of the angular structure in Mie 
theory, with the only remnant being the first side-peak in the diffraction region. The difference be- 
tween the scattering for two orthogonal incident polarizations is also much less than in Mie theory, 
and sometimes vanishes entirely. 

Other sources of experimental information include Berry (1962), Donn and Powell (1963), Kirmaci 
and Ward (1979), Sassen and Liou (1979), and the book edited by Schuerman ( 1980 ). 

Summary 

This is but a limited survey of nonspherical scattering experiments, and one oriented toward at- 
mospheric science, but it should be sufficient to give a flavor of the results. 

All the measurements to date show the following tendencies, relative to Mie theory: 

(i) damping of the oscillations vs. angle and vs. size parameter 

(ii) more side-scattering (60-120°), for unpolarized incident light 

(iii) less backscattering (140-180°) unless the particle is opaque 

(iv) nearly equal forward scattering (0-50°) 

(v) rapidly worsening agreement as size parameter increases past x = 3 to 5 

(vi) less difference between scattering resulting from the two possible orthogonal polarizations of 
the incident radiation 

A.3 THEORIES - EXACT OR NUMERICAL 

In the past 25 years, numerical techniques for scattering and absorption by variously-shaped objects 
have proliferated like weeds. This is due to the relatively few exact solutions of Maxwell’s Equations 
for the ‘simplest-case’ scattering problem (an incident plane wave). The exact solution for a circular 
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cylinder in perpendicular incidence was given by Rayleigh in 1888, and extended to oblique incidence 
by Wait in 1955 (see also Liou, 1972; Cohen and Alpert, 1979, 1980). Lorenz in the 1890’s, and in- 
dependently Mie in 1908, and Debye in 1909, gave the solution for a sphere. Yeh (1965) gave a partial 
solution for elliptic cylinders. Asano and Yamamoto (1975), in a much-acclaimed paper, gave the 
solution for a spheroid. 

It is unlikely that this arsenal will be enlarged much in the future. The reason is simple: the solution 
even for a spheroid is already so complex that it offers very little advantage over a good numerical 
s olution (except to furnish an independent check). The spheroidal solution suffers from numerical ill- 
conditioning, preventing its use above size parameter x = 30 or so (Asano, 1979); and in random 
orientation, it gobbles up monstrous amounts of computer time (Asano, 1980). Thus, even for this 
simplest of all nonspherical shapes, the exact solution behaves in many ways like a numerical solution. 

The better numerical solutions, in turn, behave in many ways like an ‘exact’ solution — for example, 
by expanding the solution in a set of orthogonal eigenfunctions. This is all symptomatic, we believe, of 
a general crumbling of the formerly rigid barriers between ‘exact’ and ‘numerical’ solutions. 

For those who wish to pursue the question of ‘exact’ solutions in more detail, we offer a few more 
remarks. The problem boils down to solving the “vector Helmholtz equation” (Morse and Feshbach, 
1953, Ch. 13) for the electric field. This is done by expanding the electric field in an infinite series of 
eigenfunctions in one of the 11 coordinate systems in which the vector Helmholtz equation is 
separable. (These series are double series in general; only for the sphere and circular cylinder do they 
degenerate to single series.) The eigenfunctions used inside the scatterer are different from those used 
outside of it, because the ‘inside’ expansion must be finite at the origin, while the ‘outside’ expansion 
must produce outgoing waves at infinity. 

Then, the usual Maxwell boundary condition (continuity of the tangential component of the electric 
field) is applied on the surface of the scatterer; the inside and outside expansions must be ‘matched up’ 
there. This is where all the difficulties arise. In every case but the few which have already been solved, 
the boundary condition leads to a set of equations which cannot be solved analytically. (If boundary 
conditions were not the major hang-up, scattering from a cube, involving only simple Cartesian coor- 
dinates, would have been solved long ago. Instead, it remains an unsolved problem.) 

Those who wish to see exact solutions for metallic scatterers in various unusual coordinate systems 
may consult the book by Bowman, et. al. (1969), where the main thrust is antenna design. Kouyoum- 
jian (1965) gives the original references for solutions for metallic cones, disks, strips, and parabolic 
cylinders. Of greater interest, perhaps, is the solution by Borghese, et. al. (1979) for an arbitrary 
cluster of spheres; indeed, it seems to this author that sphere clusters would be excellent archetypes for 
general nonspherical particles. Apparently, however, the solution of Borghese et. al. is nightmarish to 
put into practice. 


Returning now to numerical solutions, we observe that they are primarily used in the range 1 to 10 of 
size parameter (the so-called ‘resonance region’). Below 1, small-particle approximations are useful; 
above 10, large-particle approximations. Very few numerical methods have been applied for size 
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parameters above 15 or so, partly because, in general, the computational demands escalate as some 
high power of size parameter. 

Fortunately, the seemingly innumerable numerical techniques fall into two rather broad classes: dif- 
ferential equation (including finite element) methods, and integral equation methods. Each ‘new’ 
technique is usually only a minor variation on a familiar theme. It is these themes which we shall at- 
tempt to highlight here. But first we shall deal with a method that defies classification. 

Purcell- Pennypacker: aggregates of dipoles 

Purcell and Pennypacker (1977) suggested replacing a scatterer by a cubic array of point dipoles 
spaced no farther than (wavelength/4 tt) apart. Each dipole has a polarizability such that the correct 
bulk refractive index is predicted by the Clausius-Mosotti relation. Exact account is taken of all 
mutual dipole-dipole interactions, leading to a set of 6N linear equations for N dipoles. 

In this method, the surface disappears , and with it all the nastiness associated with surface boundary 
conditions. That is the beauty of the method. It is elegant in its simplicity, and therefore has attracted 
a good deal of attention. Kattawar and Humphreys (1980), for example, have applied the method to 

two nearby spheres. 

Purcell and Pennypacker themselves considered no more than 100 dipoles. Yung (1978), however, has 
been able to reformulate the method as a variational principle, enabling him to use the prodigious 
number of 16,000 dipoles (corresponding to a size parameter of about x = 6). Until computer memory 
and speed become very much larger, x = 10 is probably a practical upper limit to this method. 

Rayleigh hypothesis 

Rayleigh was interested in the scattering of a plane wave from a sinusoidal diffraction grating. Be- 
tween the peaks and valleys of the grating, there are both upgoing and downgoing waives. Rayleigh 
made the hypothesis, now associated with his name, that only the upgoing wave eigenfunctions need 
be considered in the solution. This hypothesis was accepted without question until about 1950, when it 
became embroiled in controversy (Millar, 1973). 

As it applies to scattering by a finite object, this controversy still raged through the 1960’s (Bates, 
1969; Bates, et. al., 1973; Millar, 1969; Millar and Bates, 1970). Consider a scattering object C with its 
inscribed (S in ) and circumscribed (S out ) spheres, as in Fig. A.l. At any point P outside C, but inside 
S out , there are both incoming and outgoing waves; after all, part of C is outside the sphere passing 
through P, and thus can radiate inwards towards P . . . . 

If we think in terms of rays, of course, this could only happen if C were concave. But, mathematically 
speaking, the problem exists for all objects, and the solution ought really to be of the form (see Fig. 
A.l for region designations): 


(Region I) expansion in eigenfunctions regular at the origin 
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(Region II) expansion in incoming and outgoing eigenfunctions 
(Region III) expansion in outgoing eigenfunctions. 

For some very simple problems, including Rayleigh’s sinusoidal grating, Millar (1973) has actually suc- 
ceeded in showing the quantitative limits of validity of the Rayleigh hypothesis. If the grating is 

described by 


f(z) = a cos kz 

where ‘z’ is height, then the Rayleigh hypothesis fails if ka > 0.448, or, loosely speaking, if the waves 
on the surface become too ‘violent’. By naively comparing this case to the case of Chebyshev particles 
(Eq. I of the text), and analogizing z with theta, one might expect the Rayleigh hypothesis to fail also 
for Chebyshev particles such than n e > 0.05 or so. This is curiously close to the actual limits found on 
EBCM convergence for size parameters exceeding 5 or so. 

Millar shows that a necessary and sufficient condition for the Rayleigh hypothesis to be valid is that 
the singularities of the Region III expansion all lie in Region I. Thus the hypothesis may be valid with 
one origin of coordinates, and not with another. Its validity may also change, depending on which 
eigenfunctions are used in the Region III expansion. And if the scatterer has corners or edges, this 
almost always causes singularities to occur in Region II, invalidating the Rayleigh hypothesis. 



Figure A.1. Cross-section of a general scattering 

object with surface C, inscribed sphere 
S in , and circumscribed sphere S out . 
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But following this gloomy prognosis, Millar offers a way out which seems almost too good to be true. 
We merely have to abandon our craving for an ‘exact’ solution, and satisfy the boundary conditions 
only in the least-squares sense. Then an expansion in any complete set of outgoing eigenfunctions is 
satisfactory in Region II. It will even converge uniformly in Region III (although the convergence may 
be slow, if the eigenfunctions chosen are ill-suited to the geometry). 

In the various debates over the Rayleigh hypothesis in the literature, the EBCM was always regarded 
as immune to these difficulties — apparently, because it sidesteps the whole issue of making an expan- 
sion of the scattered field in Region II. One is left with an uneasy feeling, however, that the problem 
has perhaps been papered over rather than resolved. 

Differential equation approaches 

According to Mei, et. al. (1978): 

“During the last 15 years, computation in electromagnetic scattering has been actively pursued 
almost entirely in terms of integral equations. . . .this drift to integral equations is very natural — 
and for good reasons. In integral equations, the computations are limited to the scatter er itself; 
while in the finite methods they are generally spread over the entire space. In integral equations, 
the radiation conditions are automatically satisfied, while in the finite methods they require 
special numerical treatments which are often unsatisfactory. But recent advances . . . have 
urgently demanded the results of scattering by . . . inhomogeneous bodies. The only practical ap- 
proach to such problems appears to be direct solution of partial differential equations rather 
than solution by integral equations, the formulation of which in an inhomogeneous medium is a 
difficult task.” 

This explains why methods which proceed directly from the vector Helmholtz equation, without re- 
formulating it as an integral equation, have recently regained some popularity. 

Differential-equation methods are invariably simpler in concept, and simpler in execution, than in- 
tegral-equation methods, and they avoid the singular-kernel problems in those approaches. They 
tend, however, to consume much more computer time. A well-thought-out-example is the method of 
Patwari and Davies (1966). They assume an expansion in outgoing-wave eigenfunctions in Region III. 
Then they carry this solution from S out to C (the scatterer’s surface) using a finite-difference form of 
the vector Helmholtz equation. The expansion coefficients are then determined by boundary condi- 
tions at C. 

Reilly’s (1973) method is an interesting variation. He uses a Galerkin method, and applies boundary 
conditions on S out rather than on C. 

Perhaps the simplest method of all is ‘point-matching’, reviewed by Richmond (1965). The scattered 
field is expanded in an N-term series of outgoing-wave eigenfunctions, then the boundary conditions 
are enforced at N points on C, leading to a set of N linear equations for the expansion coefficients. 
Point-matching has been used by Greenberg, et. al. (1965), among others, to study nonspherica! scat- 
tering. However, the validity of point-matching became the subject of a heated debate (Bates, 1967, 
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1969; Bates, et. a!., 1973; Millar and Bates, 1970) connected with the Rayleigh hypothesis. The general 
feeling nowadays seems to be that it has uncertain convergence, is inaccurate, and uses too much com- 
puter time. 

There have been efforts to rescue point-matching by generalizing it. In ‘least-squares matching’, one 
enforces boundary conditions at M > > N points in a least-squares sense. In ‘spectral-component 
matching 5 (Millar, 1973), one picks a set of smooth basis functions, multiplies the boundary condi- 
tion by each of these in turn, and integrates over C. (If the basis functions are delta-functions, this 
reduces to point-matching.) This can make the Rayleigh hypothesis valid even when it fails for point- 
matching. Of course, computation is increased considerably, because many surface integrals must be 
done. Spectral-component and least-squares matching are closely related from a mathematical point 
of view. 

As Mei, et. al. (1978) have noted in the quote above, the necessity to calculate scattering from in- 
homogeneous objects has breathed new life into differential equation methods. But they have revived 
mostly in the form of ‘finite-element’ methods (e.g. Morgan and Mei, 1979; Yeh and Mei, 1980), 
which Mei et. al. review in detail. 

Jntegrai equation methods 

Integral equation methods reformulate the vector Helmholtz equation as an integral equation using 
Green’s functions (Morse and Feshbach, 1953, p. 1769 ff.). The boundary conditions are included 
automatically. Within and/or on the surface of the scatterer C, one solves either for 

® the induced current J, or 

® the electric field E in . 

The scattered field is then found from the vector form of Huyghens’ principle or its equivalent. 
There are three variants, involving: 

(A) volume and surface integrals for either J or E in in one vector equation; 

(B) surface integral equations only 

(1) coupled equations for J and its spatial derivatives, with singular kernels; 

(2) coupled equations for J alone, with a higher-order singularity than in (1). 

The singularities in (B) come from the so-called ‘free-space Green’s function’. They are the cross that 
integral equation methods are forced to bear, in return for escaping the boundary condition 
nightmares of differential equation methods. 
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Livesay and Chen (1974) do nothing fancy at all. They just calculate E in by directly evaluating the 
volume integrals in a type (A) formulation. A second quadrature over E in then yields the scattered 
field. They can end up with huge matrices compared to the EBCM, but they do avoid the surface in- 
tegrals required in the EBCM. 

Some methods avoid numerical evaluation of the singular integrals in type (Bl) formulations by 
analytic methods. The EBCM (Waterman, 1965, 1971; Barber and Yeh, 1975; Mugnai and 
Wiscombe, 1980) uses different spherical harmonic expansions inside and outside the sphere passing 
through the singularity. Chu and Weil (1976) and Holt, et. al. (1978) evaluate the integrals analytically 
just in the neighborhood of the singularity. Actually, since the (Bl) singularities are integrable in the 
principal value sense, sophisticated modern quadrature techniques could probably handle them; but 
no one has tried this. 

The singularities in type (B2) formulations are non-integrable. The only approach which has so far 
been used is the Hadamard finite part idea (Senior and Weil, 1977), although the physical interpreta- 
tion of this procedure is elusive. 

“Moment methods” are a particular variant of integral equation methods (cf. the book by Har- 
rington, 1968, and the review article by Miller and Poggio, 1978). If the integral equation to be solved 
is written schematically as 


LF = G 

where L is a linear operator, then the idea is to expand F in a set of basis functions f n 

N 

F = L Mn 
n= 1 


and then to take M weighted moments of the resulting equation: 

N 

E>,( », L(f„) = f W.G 

n= 1 


Taking M < N, one then solves this set of linear equations for the coefficients a-sub-n in the least- 
squares sense. 

Moment methods are frequently applied to ‘wire structures’, either actual (like antenna arrays), or 
simulacrams of real objects. 

In spite of the revival of differential equation methods for inhomogeneous scatterers, it is worth 
noting that Wang and Barber (1979) have extended the EBCM to such objects (see also Druger, et. al., 
1979). 
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A.4 THEORETICAL - APPROXIMATE 

Approximate methods fall into the following general categories: 

• Rayleigh or Rayleigh-Gans limit 

• geometrical optics and extensions thereof 

• perturbation theory 

• thin in one dimension 

• semi-empirical 

• replace by a simpler shape 

Small particles 

Stevenson (1953) gives the general theory for scattering by arbitrarily-shaped bodies in the Rayleigh 
limit where size divided by wavelength is small. Kleinman (1967, 1978) made important improvements 
in this theory (see also Herrick and Senior, 1977; Senior, 1976, 1980). A completely analytic solution 
still eludes us, except for 5 simple shapes, including ellipsoids. The ellipsoid solution has been widely 
exploited (e.g., Huffman and Bohren, 1980). For general shapes, one must solve a simple integral 
equation for the polarizability tensor by numerical methods. 

The Rayleigh-Gans approximation is even more restricted, requiring not only small particles but a 
refractive index near unity as well (van de Hulst, 1981). The trade-off is that great analytical progress 
can be made in this case, and, like many asymptotic approximations, it is useful outside its strictly 
defined range of validity. Barber and Wang (1978) found it reasonably accurate up to refractive in- 
dices of 1 . 1 and size parameters of 1 . And there are more applications than one might imagine; for ex- 
ample. to scattering by cells in aqueous solution. 

Acquista (1976, 1980) developed an extended Rayleigh-Gans approximation, based on an iterative 
solution to Shifrin’s integro-differential formulation of scattering. It is valid all the way up to size 

parameters of 5 or so. 

Large particles 

The class of methods for particles large compared to the wavelength is extensive. Kouyoumjian (1965) 
gives an excellent review. He notes that progress in developing large-particle asymptotic solutions to 
Maxwell’s Equations (or to particular solutions like Mie theory) has been excruciatingly slow, due to 
the great mathematical difficulties (see also Kline, 1962). 

“Geometric optics” is the simplest large-particle approximation. It refers to the calculation of intensi- 
ty along ray paths as those rays experience reflection and refraction at the boundaries of a scatterer, 
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and absorption within it. The Fraunhofer diffraction pattern is normally added to the geometrical op- 
tics approximation, to account for the strong forward peak in scattering from large particles. Phase 
and polarization information can be carried along with the intensity, as van de Hulst (1981) shows, but 
usually the interference between rays is ignored, since even a slight spread in size or wavelength will 
wash it out. 

Geometric optics is particularly simple for a sphere (Liou and Hansen, 1971), because the ray paths re- 
main in a plane. For any other shape, the ray paths are twisted 3-D curves, and an analytical solution 
borders on the impossible. Instead, a Monte Carlo approach is normally taken: the scatterer is bom- 
barded with enough rays to give decent statistics on the scattered energy at each angle, lacobowitz 
(1971), Wendling et. al. (1979), and Coleman and Liou (1981) all did this for a hexagonal ice cylinder. 
Takano and Tanaka (1980) did it for circular cylinders. 

Concave particles greatly complexify geometric optics. A ray emerging from the scatterer cannot be 
logged and forgotten; this can only be done when it finally leaves the circumscribing sphere, since 
there is always the possibility that it may re-enter the particle. Thus, the ray paths may ramify almost 
endlessly. 

Reflection from a large, randomly-oriented, convex particle is rigorously identical to the reflection 
from a large sphere (van de Hulst, 1981; Hodkinson, 1963). The diffraction is close to that from an 
equal-projected-area sphere. Hence, if the particle absorbs enough to extinguish transmitted rays, 
reflection plus diffraction from an equal-projected-area sphere is an excellent approximation to the 
nonspherical scattering. 

Some authors (e.g. Pollack and Cuzzi, 1980) postulate that the transmitted rays suffer greater devia- 
tions for any non-sphere than for a sphere, accounting for the higher side-scattering from non- 
spheres. 

Chylek (1977), drawing on a result of Vouk(1948), notes that the extinction cross-section of a large, 
randomly-oriented particle always exceeds that of an equal-volume sphere. 

Maxwell’s Equations can in principle be expanded as wavelength tends toward zero (the so-called 
‘Luneberg-Kline expansion’), leading in zeroth order to the famous ‘eikonal equation’ (Kline, 1962). 
But the eikonal equation is non-linear, and the equations for the higher-order terms of the Luneberg- 
Kline expansion become increasingly complicated and non-linear; furthermore, these equations break 
down near light-shadow boundaries, where diffraction occurs. Little progress has been made on these 
higher-order equations. 

Keller (1962), in a famous piece of work, developed a theory which extends geometric optics by adding 
“diffracted rays”. These are produced whenever an incident ray hits an edge, vertex, or shadow 
boundary, and behave according to their own laws, different from ordinary rays. 
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Perturbation theory 

The idea of perturbation theory is simple enough: assume the scatterer is bounded by the surface 

r s = r 0 [1 + e f (©,</>)] 

where e < < 1. Then expand everything in powers of e. Yeh (1964, 1965) worked out the order e cor- 
rection terms for deviations from a cylinder as well as a sphere. Erma (1968a, b) corrected Yeh’s for- 
mulation and extended it to arbitrary orders. Neither author gave any numerical results. The formulas 
for scattering quantities are complicated; double series for order e, triple series for order e 2 , and so on. 

Chylek, et. al. (1978) made some first-order perturbation theory calculations and generally found 
unacceptably large errors for e > 0.10. For e = 0.05, they found roughly 1% errors compared to the 
exact spheroidal results of Asano and Yamamoto (1975). They also pointed out that perturbation 
series become useless, and may in fact diverge, once the Mie coefficients a n and b n develop sharp 
spikes (beyond size parameter x = 6 or so). This is because perturbation theory involves derivatives of 
a n and b n . When these derivatives become large, the omitted terms in the perturbation series become 
larger than the ones which are kept. 

Two particles are ‘conjugate’ if they are described by 

r s = r 0 ± g (©,</>) 

(with g such that r s > 0). Using perturbation theory, Chylek, et. al. (1979) showed that, for small g, 
the two conjugate particles have an average scattering equal to that of a sphere of radius r 0 . This is an 
interesting twist on the idea of replacing a particle by an equivalent sphere; here two particles (or a 
single “self-conjugate” particle) can be replaced by an equivalent sphere. 

Thin particles 

If a particle is thin in one of its dimensions, it is possible to approximate the integral equations describ- 
ing scattering. The work of Chu and Weil (1976), Weil and Chu (1976, 1980), and Senior and Weil 
(1977) is an example of this approach. They have applied it to both ice crystal plates and aerosol 

particles. 

Semi-empirical 

The word ‘semi-empirical’ denotes a class of methods containing a mixture of experimental and 

theoretical results. 

Emslie and Aronson (1973) and Aronson and Emslie (1980) developed such a method for the dust 
blanketing the Moon and Mars. They modeled the small dust grains as ellipsoids in the small-particle 
limit and the large grains as spheres in the geometric optics limit. Roughness on the spheres was 
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mimicked by distributing absorbing dipoles over their surfaces, which could interact with the 
geometric optics rays. Intermediate-sized particles were treated by interpolation between these two ex- 
tremes. There were probably too many tuneable parameters in this model, but it did yield agreement 
with measured emission and reflection spectra which had not been previously explained. 

An abortive semi-empirical theory which involved truncating the sharp spikes in the Mie coefficients 
a,, and b n was proposed by Chylek, et. al. (1976). The idea was that these spikes were associated with 
‘surface waves’, which would be suppressed on nonspherical particles. Unfortunately, the spikes do 
not correspond to surface waves, and in any case exact calculations by Mugnai and Wiscombe (1980), 
among others, show that surface waves persist even on rather nonspherical particles. This theory 
became embroiled in controversy (Acquista, 1978; Chylek and Pinnick, 1979) from which it has never 
recovered. 

Pollack and Cuzzi (1980) suggested a theory based primarily on the measurements of Zerull (1976). 
They use equal-volume Mie theory for nonspherical particles with size parameter x < where Xq is 
tuned in the range 3 to 10. For x > Xq, the absorption cross-section is still gotten from Mie theory, 
while the phase function is gotten from a sum of 

® Fraunhofer diffraction 

• reflected rays from a sphere 

® transmitted rays fitted to mimic Zerull’ s measurements (involving a second tuneable 
parameter). 

The single-scattering co-albedo of the sphere is multiplied by the ratio 

surface area of particle/surface area of sphere 
to get the co-albedo for the particle. 

Coletti (1984) has devised another semi -empirical theory based on his own measurements, and similar 
in some respects to that of Pollack and Cuzzi. 

Replace by a simpler shape 

Perhaps the simplest approximation is to replace the particle by a simpler shape, which however has to 
fulfill two sometimes contradictory requirements: 

® the scattering properties of the simpler shape have to be relatively easy to calculate, and 

• they have to resemble those of the actual particle. 
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Spheres are by far and away the most popular replacement. The only question then is: of what 
radius? Equal-volume spheres are the choice of 90% of all investigators. Equal-surface-area 

ana equal-projected-area spheres also have their devotees. A less well-explored option is to match, 
not a material property of the particle, but an optical property: extinction cross-section, for example. 

The argument for equal-volume spheres (aside from the relative ease of measuring particle masses) is 
that, at least for particles of small size parameter, scattering depends primarily on volume, not shape. 
This is gainsaid, however, by Aronson and Emslie (1980) and by Huffman and Bohren (1980), who 
find that, even when one averages over shape and orientation, ellipsoids do not behave like spheres. 
Likewise, Liou (1974) has found that equal-volume spheres are a poor substitute for ice cylinders in 
the 8-12 micron infrared window region. 

Equal-projected-area spheres are generally regarded as the best replacement for particles of large size 
parameter, because the forward diffraction peak depends primarily on projected area (Hodkinson, 
1963). Since this peak is the prime determinant of the low-order phase function moments, and since 
these moments are the prime determinants of multiple scattering, it is worth considering the equal- 
projected-area sphere for that application. Holland and Gagne (1970), among others, find that the 
scattering from 0 to 40-50° is best approximated by equal-projected-area spheres. 

A.5 SPECULATIONS AND IDEAS 

Nonspherical particle scattering is so exceedingly rich and various that there is no hope of jumping 
immediately, by a sort of Aristotelian contemplation, to an all-embracing theory. The rush to 
‘semi-empirical’ theories is, in our opinion, premature; it has already led to one bad mistake, and is 
likely to lead to more. Instead, we should proceed inductively, examining many special cases first. 

Mie theory should be abandoned only grudgingly. Single equal-volume spheres are surely not an ade- 
quate replacement for all, or even most, nonspherical particles. But we are far from having exhausted 
Mie theory. There are still several largely unexplored replacement possibilities: 

® a sphere having the same volume-to-surface-area ratio 

• a size distribution of spheres (e.g., Wang et. al., 1979) 

• a multi-layer sphere (e.g., Kerker, 1969, Ch. 5) 

• a sphere with continuously-varying refractive index (e.g., Kerker, 1969, Ch. 5). 

The last two possibilities are based on the simple idea of trading off nonsphericity against in- 
homogeneity. Imagine a tumbling nonspherical particle as a ‘fuzzy sphere’, with a solid core shading 
gradually outwards. Between the core and the periphery, the refractive index would gradually change 
to that of the surrounding medium. By picking a functional form for this change that allowed a 
reasonably simple radial solution, with one or two adjustable parameters, it might be possible to 
match nonspherical scattering properties. Because scattering from an inhomogeneous sphere is orders 
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of magnitude easier to calculate than scattering from any non-sphere, this idea is immensely 
attractive. 

Modifications to Mie theory to account for non-sphericity may also be possible; but not until a 
definite physical interpretation can be assigned to the terms. In this regard, Complex Angular Momen- 
tum Theory (Nussenzweig, 1979) may prove useful for larger particles, since it replaces the Mie expan- 
sion with a series of only a few terms, each of which can be assigned a simple physical meaning. 

If we are finally forced to abandon Mie theory for some particularly recalcitrant cases, we should at 
least try to use equivalent spheroids . This gives us two parameters to tune to match the actual 
nonspherical results, and that may be enough. 

Concave particles, and particles with voids, may defy all our attempts to replace them with equivalent 
spheres or spheroids. Since such particles are of great practical importance, we should make every ef- 
fort to find their scattering properties, both experimentally and from exact theory. Our own EBCM 
calculations have shown great differences between the phase functions of convex and even mildly con- 
cave particles. 

We need to develop a minimal, prototypical set of shape parameters. Up to now, we have either pick- 
ed shapes which were easy to calculate for, or which Nature thrust upon us. Very little thought has 
been devoted to the abstract concept of ‘Shape’, and what characterizes it. Greenberg (1980) has sug- 
gested some of the shape information which may be important. Our own list would include: 

® some measure of surface roughness 

• the 3 semi-axes of the ellipsoid which best approximates the particle 

• the hole volume or ‘porosity’. 

Present numerical methods need to be speeded up dramatically. Even the fastest of them, like our vec- 
torized EBCM code, consume far too much computer time. At the same time, faster ways of averag- 
ing over orientation need to be found; this is a real roadblock in present scattering calculations. 

Once faster numerical methods are available, it will become possible to routinely average over shape, 
size, and orientation, something which is now impractical except in test cases. Then we can more in- 
cisively answer questions as to how much each kind of averaging reduces spherical-nonspherical 
differences. 
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APPENDIX B 


Parallel and perpendicular scattered intensities vs. scattering angle for two fixed orientations, nose-on 
(dotted line) and perpendicular to nose-on (dashed line); for Chebyshev particles T 3 ( + 0.10), 
T 2 ( — 0.10), T 4 ( + 0.10), and T 6 ( — 0.10); and for size parameters x = 5 and 10. Spherical results (solid 
line) are shown for comparison. Note the much larger sphericai-nonspherical differences in fixed 
orientation than in random orientation. These results are intended for direct comparison with 
experiment. 


142 


SINGLE SCATTERING FROM NONSPHERICAL CHEBYSHEV PARTICLES 



A.USN31NI nanvyvd 


Appendix B. Nonspherical SCATTERED INTENSITIES (parallel and perpendicular) vs. 
ANGLE for Chebyshev particles in two fixed orientations 




lel=0. 10 x=5. t? = 0.,90. T 3 + lel-0.10 x=5. t?=0.,90. 
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Phase functions in the side- and back-scattering regime (60-180 degrees) vs. scattering angle for ran- 
domly oriented Chebyshev particles T 3 , t 2 , t 4 , t 6 , and T 8 ; for various deformation parameters be- 
tween - 0.20 and 0.20 (as shown above each plot); and for size parameters x = 2, 3, 4, 5, 6, 8, 10, 15, and 
20 . Solid line is for Ax = O.lx size-averaged spherical result, interrupted lines are for non-sheres: dot- 
ted, dashed, and dot-dash lines refer to the various e values shown, in increasing order (in genera!, the 
larger the value of | e | , the greater the general deviation from the spherical result). Not all particles are 
represented at larger values of x because of an inability to achieve satisfactory convergence of the 
EBCM. 

Next to each phase function plot is a plot of the percentage deviations of the nonspherical phase func- 
tions from the spherical one, but for the full range of angles (0-180 degrees) rather than just 60-180 
degrees. The interrupted lines in these plots refer to the same nonspherical particles as in the adjacent 
phase function plot. Note that when the percentage differences become larger than about 100% the 
vertical scale becomes logarithmic, and in order to display both positive and negative values, percent- 
ages smaller than 1% in magnitude are assumed to be zero. 
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Degree of polarization vs. scattering angle for randomly oriented Chebyshev particles t 3s t 2 , t 4 , t 6 , 
and T 8 ; for various deformation parameters e between -0.20 and 0.20 (as shown above each plot); 
and for size parameters x = 2,3,4,5,6,8,10,15, and 20. Solid line is for the Ax = O.lx size-averaged 
spherical result, interrupted lines are for non-sheres: dotted, dashed, and dot-dash lines refer to the 
various e values shown, in increasing order (in general, the larger the value of | e | , the greater the 
general deviation from the spherical result). Not all particles are represented at larger values of x 
because of an inability to achieve satisfactory convergence of the EBCM. 
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<V> lel=0.05.0.10 x=2. 
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(V) l£l = 0. 05,0. 10,0. 15 X = 6. 
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(T 3 “) lel = 0.05,0. 10 x=20. 
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